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FOREWORD 


This document was prepared by the Engineering Division of the Pratt and 
Whitney Group of United Technologies Corporation East Hartford, Connecticut, 
to describe the selection, incorporation, and evaluation of the best available 
finite difference scheme to reduce numerical error in the Combustor 
Performance evaluation codes. This task was accomplished under the Mational 
Aeronautics and Space Administration (NASA) sponsored Error Reduction Program, 
Contract NAS3-23686. The work was performed under the direction of NASA 
Project Manager, Mr. Russell W. Claus. 

This report has been assigned the Contractor's report number PWA-5928-25. 
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1.0 SUMMARY 


i he National Aeronautics and Space Administration sponsored a program to 
select, incorporate and evaluate the best available finite difference scheme 
to reduce numerical error in combustor performance evaluation codes. This 
report describes the details of this study. 

The combustor performance computer programs chosen for this study were the two 
dimensional and three dimensional versions of Pratt and Whitney’s TEACH code. 

The criteria used to select schemes required that the difference equations 
mirror the properties of the governing differential equation, be more accurate 
than the current hybrid difference scheme, be stable and economical, be 
compatible with TEACH codes, use only modest amounts of additional storage, 
and be relatively simple. 

The methods of assessment used in the selection process consisted of 
examination of the difference equation, evaluation of the properties of the 
the coefficient matrix, Taylor series analysis, and performance on model 
problems. Five schemes from the literature and three schemes developed during 
the course of this study were evaluated. This initial evaluation resulted in 
the selection of the two most promising schemes. Quadratic Upwind Differencing 
Scheme (QUDS) and Bounded Skew Upwind Differencing Scheme Two (BSUDS2), for 
incorporation into 2D-TEACH for further evaluation. The accuracy and stability 
of these schemes were assessed by using laminar and turbulent flow test cases. 
These test cases, although two dimensional, contain important flow features 
found in gas turbine combustors. During the evaluation, it was found that QUDS 
was unstable and, hence, BSUDS2 was selected for incorporation into 3D-TEACH. 
This scheme was further evaluated by using a 3D-test case, modeling of a jet 
in cross flow. 

This effort resulted in the incorporation of a scheme in 3D-TEACH which is 
usually more accurate than the hybrid differencing method and never less 
accurate. It is expected that overall improvement in accuracy resulting for 
complex flows will justify the increased cost of using this scheme. However, 
this study can only be considered as a first step in the process of developing 
the most suitable scheme for combustor performance codes. A number of 
questions have been generated as a result of this work and it is expected that 
answers to these questions will lead to further improvements in the accuracy 
and stability of the BSUDS2 scheme. 



2.0 INTRODUCTION 


2.1 OBJECTIVES 

The main objective of the NASA-sponsored Error Reduction Program was to 
select, incorporate and evaluate the best available technique for the 
reduction of numerical diffusion in a 3D combustor performance evaluation 
code. The study focused on improvements in accuracy of computer programs of 
the TEACH (for Teaching Elliptic Axisymmetric Characteristics Heuristically) 
variety which were developed originally by Professor A. D. Gosman and 
co-workers of Imperial College., London, (e.g. Ref. 1 ) and are in general use. 

The need for such a program can best be appreciated by comparing the numerical 
and exact solutions for the spreading of a passive scalar (dye in water for, 
example) in a simple flow field, Fig. 2-1. The differencing scheme used for 
this computation is the hybrid method which is used in almost all TEACH 
combustor codes. It can be seen that even for this simple flow situation, 
numerical (or artificial) diffusion generated by the differencing scheme 
greatly smears the profile. Numerical diffusion, especially in 
three-dimensional versions of TEACH, can become so large as to obscure the 
effects of the turbulance model in turbulent flow calculations. The accuracy 
of these codes becomes dependent on the number and distribution of nodes used 
in the finite-difference calculations; i.e. s it is difficult to achieve a 
"grid-independent" solution. The number and distribution of grid nodes varies 
from user to user so that the results of parametric studies requiring more 
than variations in boundary conditions are questionable. This places a major 
restriction on the utilization of these codes for the design and development 
of practical combustors. Hence, there is a requirement for improving the 
accuracy of the differencing scheme presently embodied in these codes. This 
problem is addressed by the current program. 



Y, MM 

Convection and diffusion of a scalar step profile in constant property 
uniform inclined flow, cell peclect number, Pex = Pey = 60 
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Figure 2-1 Need for Error Reduction 



2.2 APPROACH 


The strategy adopted in this study was to select a number of candidate schemes 
and assess them for accuracy, solution stability, and overall 
cost-effectiveness. The accuracy assessment included consideration of the 
conservation and boundedness properties of the schemes, as well as 
conventional Taylor series error analyses and application to test cases. 
Solution stability was assessed heuristically by examining the properties of 
the coefficient matrix which each scheme generates. Cost-effectiveness was 
judged on the combined outcome of the foregoing assessments. 

Based on this initial evaluation, two schemes were selected for incorporation 
into Pratt and Whitney's two-dimensional combustor performance code, 2D-TEACH 
(This computer program is the two-dimensional version of the three-dimensional 
code, 3D-TEACH, described in detail in Appendix B). The stability, accuracy 
and cost effectiveness of these two schemes were then examined by running a 
number of laminar and turbulent flow test cases. These test cases demonstrate 
typical flow features encountered in gas turbine combustors. 

The more suitable of the two schemes was then incorporated into Pratt and 
Whitney's three-dimensional combustor performance code, 3D-TEACH (Appendix B). 
The improvement in accuracy of. this scheme was evaluated for a test case 
modeling a row of jets in a cross flow and comparing the calculations against 
both experimental data and calculations performed using the hybrid 
differencing scheme. 

2.3 ORGANIZATION 

In Section 3 of this report a brief description of the sources of errors in 
the present combustor performance codes, 2D- and 3D-TEACH, is given and 
reasons for restricting attention in the present study to errors caused by the 
differencing scheme are explained. In Section 4 the procedure for selecting 
the two most promising schemes for incorporation into 2D-TEACH is described. 
This procedure involves selection of a number of schemes from the literature 
and evaluating their cost effectiveness by assessing their accuracy, 
stability, complexity, storage requirements and compatibility with 2D- and 
3D-TEACH. In Section 5 implementation of the selected schemes into the TEACH 
codes is detailed and in Section 6 a description of the two- and 
three-dimensional test cases is given and the results of the computations are 
discussed. In Section 7, concluding remarks are given and in Section 8 
recommendations for future work are presented. 
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3.0 BACKGROUND 


An error reduction program in computational fluid mechanics faces the problem 
of selecting which of several errors are to be reduced since there are many 
sources of numerical error present in current design analysis codes. In this 
section, these error sources are described and the reasons for concentrating 
on errors due to the finite - difference method are given. 

3.1 SOURCES OF ERROR IN PRESENT DESIGN CODES 

The sources of error in present computer programs can be explained by 
referring to Figure 3-1 showing the four major steps of a typical 
computational procedure. In the following sections the errors introduced in 
each of these steps will be discussed with reference to the structure and 
models incorporated into the 2D and 3D-TEACH codes. 

3.1.1 Assembly of Equations 

In the first step, the governing partial differential equations are assembled. 
These consist of the Navier-Stokes equations for mass and momentum 
conservation, the energy equation, and the species transport equations 
together with additional auxiliary equations. The governing partial 
differential equations can be regarded as exact. The auxiliary equations, such 
as those used to represent various physical processes like chemical reaction, 
turbulence generation, etc., are generally only approximations in part because 
the relevant physical processes are known only approximately. In subsequent 
steps, other approximations (e.g,, the finite-difference representation of the 
governing equations) introduce additional errors. 

For example, the instantaneous value of a dependant variable, 0 , in a 
turbulent reacting flow is usually taken as the sum of a time-mean value and a 
fluctuating value, a computational convenience, whose physical realism in 
increasingly questioned. When this definition is introduced into the 
Navier-Stokes equations and these are then time-averaged, further 
approximations to simplify the resulting Reynolds-averaged equations are 
Introduced; specifically, fluctuations in laminar viscosity and density are 
often neglected. 

As another example, infinitely-fast reaction rates are often assumed so that 
chemical reactions are represented by the one-step, irreversible reaction: 

Fuel + Oxidant-MProducts (3.1) 

In this case, the combustion rate is assumed to be controlled by the turbulent 
mixing of eddies containing the reactants. From a chemical kinetics point of 
view, this reaction scheme is extremely crude and is relevant only in 
situations in which kinetics do not control the heat release rate. 

Fortunately, most of the gas turbine combustor operating envelop is in this 
category , 
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3.1.2 Physical Modeling 

It has been established (Ref. 1 ) that a hierarchy of physical models exists. 

This hierarchy consists of models for: 

o Turbulence 

o Fuel -Spray Vaprization and Distribution 
o Combustion 
o Thermal Radiation 



Figure 3-1 


Flow Diagram of Computat 


ional Procedure 
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Clearly, there is little advantage to paying the computational cost of using, 
say, a radiation model of higher order accuracy than any of the models 
preceding it in the hierarchy. 

T he Re ynolds averaged and other equations contain terms such as 7 ujui and 
p ujm where m is the mass fraction of chemical species "i". Modeling is 
required to provide expressions for these correlations in terms of either 
known or calculable quantities. 

Reynolds stresses can consist of two parts - a shear stress and a normal 
stress. The normal stress is obtained simply from the fluctuating dynamic 
pressures, while Boussinesq's analogy is used to relate the shear stress to 
the velocity gradient through an eddy viscosity pt. Thus variable density 
flow the Reynolds stresses can be expressed as: 


-p 


/ 3u ■ 3u- 3u. \ — - 

u'u' - IV — ± + — 1 - 2/3 — - 6. •) - 2/3 p K 6^ 

1 J C l Skj 3x k U j ‘J 


( 3 . 2 ) 


where K is the turbulent kinetic energy and is defined as: 


K = 1/2 (u' 2 + v' 2 + w' 2 ) 


The eddy viscosity, M t is obtained by dimensional arguments from the 
Prandtl- Kolmogorov definition. 


where 


- c K 

u t = p — — 


( 3 . 3 ) 


e = dissipation rate of K 


The eddy viscosity is evaluated by means of transport equations for K and e. 
using the familiar two-equation turbulence model. 

Since the flow field is calculated using an effective turbulent eddy 
viscosity, it is computationally convenient to base the turbulent heat and 
mass transfer rates on effective thermal and mass eddy diffusivities. 
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In general, the eddy diffusivity is used to calculate the flux of a scalar by: 


‘ u i 0 


— r 3 0 


(3.4) 


where. 


0 = scalar quantity such as temperature or species concentration 
It 0 = turbulent eddy diffusivity for 0 

The eddy diffusivity is found from the ratio of turbulent kinematic eddy 
viscosity v-f- to the turbulent Prandlt or Schmidt number, c^, for 0 . 

The turbulence model is a second order, mean-field closure to the equations. 
For the eddy viscosity approach, the two-equation model is the most general 
and sophisticated representation, and it is not computationally expensive. 
Sophistication comes from the use of differential equations to describe both 
the velocity scale and length scale to which the eddy viscosity can be 
related, rather than relying upon an a priori length scale specification as 
used in the mixing-length approach. 

There are a number of limitations to the two-equations turbulence model. It is 
well known that the assumptions used to derive the turbulence kinetic energy 
dissipation rate (e) equation are somewhat arbitrary. More importantly, use of 
the gradient diffusion idea itself has long been challenged. There are 
objections to the assumption that the Reynolds stresses depend on only the 
local mean rate of strain as well as to the assumption that the stresses are 
proportional to the local rate of strain. The "constant" of proportional ity 
really depends on the ratio of local production and dissipation of turbulence 
energy, and this ratio is not actually a constant. A further weakness is the 
adoption of a single velocity scale at a point in the flow although it is 
known that this scale can vary from point to point. The implication of a 
single scale is that the turbulence is isotropic. Turbulent flows usually 
posses some degree of anisotropy and some flows (e.g. flows with swirl or with 
large streamline curvature in the streamwise direction) produce turbulence 
that is highly anisotropic. The velocity and length scales have to be the same 
order of magnitude as the mean field motion. This is only true for flows 
dominated by simple shear forces; buoyancy forces, for example, have separate 
scales. It is implied that the turbulent motions have a small scale compared 
to that over which the concentration of a diffusing quantity changes 
significantly, yet most of the larger eddies in a turbulent flow do not 
satisfy this condition whether they are coherent or not. Thus, material can be 
transported by vortical motion against the gradient of the scalar. Williams 
and Libby (Ref. 2) have called this process "counter-gradient diffusion," 
while Spalding has used the more descriptive phrase "pressure-gradient 
diffusion" (Ref, 3). By relying on local mean rates of strain consideration 
of the effect of flow history on turbulence structure is lost. 
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It can also be appreciated that there are levels of approximation that are, 
introduced into the averaging process that results in an enormous 
simplification for turbulence modeling. Similar approximations must be used 
for all the physical modeling used. 

3.1.3 Computer Solution 

When both the governing and auxiliary equations are assembled, they form a 
simultaneous set of non-linear partial differential equations and algebraic 
equations. 

Numerical solution of the equation set is necessary. Conventional numerical 
methods available to solve equation sets of these types can be broadly divided 
into finite difference and finite element methods, although the dividing line 
is not distinct. Finite difference methods have a considerable background in 
the fluid dynamics area, and most solution approaches, including TEACH, 
utilize finite differences. 

The finite difference analog of the differential equations is obtained by 
overlaying a computational mesh on the flow domain, and obtaining the basic 
finite difference form of the partial derivatives for every node of the mesh 
from a control volume approach (Ref. 4 ). The finite difference expressions, 
when substituted back into the differential equations, yield a set of 
linearized, algebraic equations for every node of the mesh. Thus, there are as 
many sets of equations as there are nodes in the calculation domain. These 
sets, along with the boundary conditions for the problem, can then be solved 
to give solutions for the entire flowfield. 

The accuracy of a differencing scheme can be judged from the qrder of the 
terms of an equivalent Taylor Series that have been retained in the expansion. 
Unfortunately, the requirements of numerical stability are opposite to those 
of accuracy with respect to these terms. Achieving a balance between accuracy 
and stability can be particularly trying in the case of a chemically reacting 
flow because of the coupled nonlinearities which exist between the chemical 
and fluid mechanical processes. The spatial differencing of the convective 
terms of the conservation equations in an Eulerian coordinate system can 
result in numerical diffusion. Use of a higher order differencing scheme 
eliminates or significantly reduces this diffusion. However, the use of 
central -differencing method, for example, often produces oscillations in the 
solution that have no physical significance (Ref 5 ). The use of an upwind 
differencing, or donor cell, technique eliminates oscillations but introduces 
a diffusion-like term into the difference equations. Thus, while "numerical 
damping" suppresses oscillation, it leads to significant additional diffusion 
of the convected parameter. For flows with combustion, these parameters might 
be species concentration, temperature, etc. Unfortunately, diffusion of these 
quantities is responsible in a physical sense for flame propagation. 

Therefore, a severe restriction can be placed on the quality of quantitative 
prediction (Ref. 6 )• 
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It can be argued that use of upwind differencing in regions where convection 
strongly dominates streamwise diffusion is reasonable since the local upstream 
values of the field variables are swept downstream virtually unchanged, 
whereas in high-diffusion regions the form of the relatively small convection 
terms is not important. In regions where the two transport mechanisms are 
comparable, a switch to more accurate central differencing for convection or 
use of a suitably weighted combination of central and upstream differencing 
can be used. 

This somewhat narrow view of complex flows has led to the appearance and use 
of a popular and successful hybrid central /upwind differencing scheme 
(Ref. 7). This scheme is currently used in TEACH codes. The method uses 
central differencing for convection and diffusion fluxes when the absolute 
value of the Peclet number for the control volumes existing about grid nodes 
is less than, or equal to, two: upwind differencing for convection fluxes and 
neglect of diffusion fluxes is used otherwise. Peclet number defines the 
relative importance of convective and diffusive transport and is numerically 
equivalent to a cell Reynolds number. 

To use successfully the hybrid differencing scheme for complicated flows, care 
must be taken in establishing the computational grid upon which the 
calculations are performed. The approximations of the algebraic expressions 
used to represent the partial differential equations becomes asymptotically 
exact as the distance between the grid nodes is reduced. In the limit, the 
number of nodes can be increased until an asymptote to the solution of the 
differential equations is achieved. In practice, this increase is limited by 
computer storage and the cost of the calculation. However, it is not only the 
number of nodes that are used in determining the accuracy of a solution which 
is important, but also the distribution of those nodes within the flowfield to 
be determined (Refs. 8 , 9 ). This nodal distribution is important because 
whenever curvature of the flow in the streamwise coordinate direction exists, 
a truncation error arises in the solution (Ref. TO ). In addition, there is 
also a problem in multidimensional flows of streaml ine-to-grid skewness 
(Ref. 11 ). With upwind differencing, these effects start to have a damaging 
effect on solution accuracy when the Peclet number exceeds two. 

It has been concluded (Ref. 1 ) that the hybrid finite differencing scheme, 

although yielding physically realistic solutions in all circumstances, 
introduces excessive numerical diffusion for many two-dimensional flows and 
for all three-dimensional flows because presently available computer storage 
is generally not sufficient to permit local adjustment of the grids as 
described above. Thus, solution accuracy is presently controlled by the 
numerics rather than the hierarchy of physical modeling. 

3.1.4 Representation of Geometry 

A drawback of the present calculation methods based on TEACH- type computer 
programs is the lack of flexibility with respect to irregularly-shaped 
boundaries for the calculation domain. Therefore, the geometry is 
"discretized" to fit the coordinate system. 
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Figure 3-2 shows a typical modern^ annular combustion chamber for an aircraft 
gas turbine engine. To calculate flows in such a combustor using a TEACH code 
meahs that curvilinear surfaces must be represented using "stair-steps." 
Figure 3-3 gives a two-dimensional (axi symmetric) example of such a 
representation. 





Figure 3-2 Cross Section of a Typical Modern Annular Combustion Chamber 
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The use of stair-step geometries has a number of implications: First, surface 
areas are not correct. Thus, irrespective of physical modeling and numerical 
accuracy, calculation of wall shear stress and surface heat transfer rates can 
never be correct. Second, adequate representation of the geometry bounding the 
flow to be calculated usually requires more computer storage than is available 
on the present generation of computers. Mesh refining to control numerical 
diffusion is not therefore possible, and the calculated flowfield may be 
influenced incorrectly by the geometric representation. 

3.1.5 Solution Algorithm 

By their very nature, computational fluid dynamies (CFD) codes are large 
consumers of computer resources. The cost of a CFD calculation depends on the 
number of grid nodes, the number of equations to be solved at each node, and 
the rate of convergence produced by the solution algorithm. Generally, the 
number of nodes and equations are determined by the scope of the problem under 
investigation. If the convergence rate is too low, considerations of cost and 
economy may result in termination of the solution, even though the residuals 
in the relaxation process might still be high. An arbitrary limit is thus 
placed on solution accuracy by considerations of cost. 

Low convergence rates in existing codes can occur, typically, if the 
calculation domain is large and has several entering streams so that the 
flowfield thereby contains a number of recirculation zones. In addition, an 
increased number of grid nodes is necessary to resolve adequately all the flow 
features. For such strongly elliptic flowfields, the weak coupling of the 
momentum and continuity equations through the SIMPLE (Semi -Implicit Method for 
Pressure-linked Equations) algorithm and the weakly implicit nature of the ADI 
(Alternating Direction Implicit) matrix solution procedure can result in low 
rates of convergence. In addition, such a flowfield can become physically 
marginally stable, nonstationary, or bistable in character, producing 
eddy-shedding, Coanda effects, separation, etc. A combination of entering 
flows could exist for which the presumed steady-state solution does not exist. 
The flow pattern may then change from iteration to iteration and the 
convergence rate may become unacceptably low. For nonstationary flows, 
divergence can result. 

The equations to be solved can also present a convergence problem as, for 
example in reacting flows with the chemistry modeled using "stiff" exothermic 
reaction rate expressions. Similar problems arise in the turbulence equations 
(Ref, 12) when the equations are solved iteratively and uncoupled. Under 
such circumstances, severe under- relaxation is usually required to achieve 
convergence and the time to achieve a solution becomes unacceptably high. 

The finite differencing used to approximate the partial differential equations 
generates the coefficient matrix for the equations. Matrix conditioning can 
influence convergence. If the coefficient matrix is always diagonally 
dominant, then any fast matrix solver may be used without difficulty. However, 
if diagonal dominance is lost, then the ADI method tends to be unstable; if 
the difference method produces negative coefficients, divergence usually 
results. Therefore, the solution accuracy depends not only on the order of the 
differencing scheme, but also on its compatabil ity with the solution algorithm. 
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3.2 UTILITY OF PRESENT DESIGN CODES 

From the above discussion, it appears that the numerical inaccuracies of 
present design codes limit their usefulness to the designer. 

However, these codes are adequate for many engineering applications. The only 
practical engineering alernative is to conduct numerous experiments verifying 
a design. With the cost of fuel, materials and manpower increasing rapidly, 
this alternative is becoming an increasingly expensive proposition. Hence, 
even an approximate answer that eliminates some testing is acceptable. 

Although these codes lack quantitative accuracy, their qualitative accuracy 
has been demonstrated in a number of test cases. This qualitative accuracy 
allows an engineer to conduct parametric studies with confidence, allowing 
quick preliminary screening of design ideas. This ability is also useful in 
diagnosing and solving development problems. Once the problem is simulated 
with the code, it is then usually much quicker to develop an acceptable 
solution by using the computer than a rig. 

3.3 REASONS FOR CONCENTRATING ON THE DIFFERENCING SCHEME FOR ACCURACY 
IMPROVEMENT 

Several sources of error in the present design codes were described in Section 
3.1. Although these codes are still useful for design purposes, improvement in 
quantitative accuracy will increase their utility. The question then has to be 
asked: how can the accuracy of these codes be improved. One strategy is to 
develop a new generation of codes that eliminate the statistical description 
of turbulence, introduce subgrid turbulence scaling and combustion models, 
reduce numerical diffusion, use a body fitted coordinate system, and use a 
faster solver than the ADI scheme. This strategy requires a long lead time 
because most of the models are in the development stage. In addition, it takes 
several years to turn a research code into a production code. This approach is 
best suited for a university and several universities are already working on 
different aspects of this strategy. 

Another strategy to improve the present production codes in a relatively short 
period of time is to work on only one aspect of the problem. This strategy can 
yield only a limited improvement in accuracy when compared to the potential 
benefits of developing a new generation of combustor design codes; however, if 
only mature models are incorporated into the calculation procedure, the 
chances of success are large. It was shown in Section 3.1 that solution 
accuracy is presently controlled by numerics rather than the hierarchy of 
physical modeling. Hence, improvement in the differencing scheme is the first 
area that needs to be investigated and, if successful, will have an immediate 
impact on the accuracy of the code. There are several differencing schemes 
available in the literature that have been tested by a number of investigators 
and have shown promise in model problem studies of some simple flows. It is 
thought that these schemes have reached a level of maturity that they can now 
be tested in modeling more realistic flows. It was with these considerations 
in mind that NASA sponsored the Error Reduction Program to select, 
incorporate, and evaluate an improved-accuracy finite difference scheme in 
3D-TEACH. 
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4.0 SELECTION OF MORE ACCURATE FINITE-DIFFERENCING SCHEMES 

Having established in the previous section that improvements to the present 
finite-differencing scheme can yield significant benefits in the short term, 
the task then is to develop such improvements. During the development process, 
a number of constraints must be kept in mind: 

1. The scheme is to be implemented in TEACH- type codes used to design gas 
turbine combustors; these codes are operated by engineers who are not 
necessarily expert in computational fluid mechanics. Hence, the scheme has 
to be robust, require no attention from the user, and yield results that 
are always physically plausible. 

2. Since the selected scheme is to replace the hybrid scheme without any 
other changes being made to the code, it is important that the selected 
scheme be compatible with the other parts of the code such as the solution 
algorithm. 

3. The scheme should be capable of computing accurately flows of the type 
that occur in . a gas turbine combustor; thus, testing of candidate schemes 
must include some model problem studies representing realistic flowfields. 

Several improved finite-difference schemes were selected from the literature 
and subjected to an initial screening process. Two schemes were then selected 
that were believed to be capable of calculating realistic gas turbine 
combustor flowfields, and were compatible with the present code. These two 
schemes were then incorporated into the 2D version of the TEACH code. After 
using the revised computer program to calculate several laminar and turbulent 
flow test cases containing important flow features common to gas turbine 
combustors, the more promising scheme was selected for use in the 3D version 
of TEACH. The initial screening process and selection of the two most 
promising schemes are described in this section. 

4.1 CRITERIA OF ASSESSMENT 

A scheme to be implemented into the present design codes should be not only 
more accurate than the present hybrid difference scheme but it should satisfy 
certain other criteria which are discussed below. 

4.1.1 Mirror Differential Equation Properties. 

Before any scheme is judged for accuracy, it is necessary to ensure that the 
approximation to the discretized equation mirrors certain key properties of 
the original differential equation. This discussion will be facilitated if a 
general discretized form of the equation to be solved is derived. The 
two-dimension form is presented, but similar remarks apply to the 
three-dimension case. 

A prototype transport equation for a general scalar entity, <t > , which may stand 
for a velocity component, temperature, concentration, turbulence energy, etc, 
can be written as: 
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(4.1) 


9X ( ' ,U!l ' r w’ ' T $ - S 

where p and r are the density and diffusivity, respectively, S is the local 
volumetric source (sink) rate and u and v are the velocity components in 
directions x and y, respectively. 

Integral forms of the above equations can be written for finite regions. If 
one integrates Equation 4.1 over the region defined by the dashed lines in 
Figure 4-1, one obtains: 



where w, e, s and n denote the four surfaces of the region. 



Figure 4-1 Region of Integration Defined by Grid Lines 
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Since the majority of the schemes examined are based on the integral analysis, 
(usually referred to as the "finite volume method" or "FVM", Ref 4 ) it will 
be useful to assemble and examine a relatively general form of discretized 
transport equation derived in this way. This analysis will follow the general 
lines of Gosman and Lai ( 1982) (Ref. 13 ). In the FVM, attention is focussed on 
the integrals of Equation 4.2, those on the left-hand side of which represent 
total transports by convection and diffusion through the cell face in 
question. If the transport through the w face is denoted by F w . then: 

^n f , 

Fw * y j (pU£i ~ r W | x w d * (4.3) 

•9 

which may be written, using the Mean Value Theorem (MVT) of the differential 
calculus, as: 


rid 

F w = [< ^>w • (r w'w ] iy (4.4) 

where the subscript w now denotes an average along the cell face and 
Ay = yn~ys 1S the face "area". A more compact and convenient expression is 
obtained through further use of the MVT to give: 


F 

w 


C d 


w w 


r Ay 

1 w J 


(—) 
3X W 


(4.5) 


where C w = (pu) w Ay is a convection coefficient; C w and r w are obtained 
by suitable averaging of the p, u, and rfields. 

It is at this stage that the major approximations are introduced, the purpose 

of which is to relate d w and ( a0/ax) w to d values at the grid intersec- 
tions or "nodes", which are labeled in the point-of-compass fashion as F, N, 

S, and W in Figure 4-1. Here it is assumed that d w is defined in terms of 

nodal values di by: 


K - i °w 

where the <2^ are weighting factors andS w denotes summation over 
specified nodes in the vicinity of w. In a similar way (3d/ ax) will be 
approximated by: wi 

(«) w * (S + + 8 i+ + ^ - 5-- B w-~ «*,■) /4x 

where the s w are further weighting factors, w + and w _ denote locations 
to the right and left of the cell face, respectively, and ax = Xn-xy. 
Insertion of the above expressions into Equation 4.5 finally yields: 


F - C 2 a 1 ' ~ d (2. b 1 d. - j b 1 d,) 

w w w i w - i 7 - , r 1 

w w+ w+ w— w- 


(4.6) 


(4.7) 


(4.8) 


where d w = r w Ay /ax 
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The discrete analogue of the integral transport Equation 4.2 is obtained using 
F w from above and analogous formulae for the remaining F's and substituting 
into the equation: 

F w~ fr e +F s -f: n + Su = 

Before this is done it will be useful for illustration purposes to specify a 
particular set of participating nodes, as this will allow the properties of 
the resulting equation to be more readily perceived. Here one chooses, for the 
w face, the "nearest- neighbor" set SW, W, NW, S, P and N (see Figure 4-1) 
which, together with analogous sets for the remaining three faces give rise to 
the "compact" nine-point molecule illustrated in the Figure 4-1. The resulting 
equation, when arranged into a substitution formula for <5 p , runs: 

Vp *•§ Vj + Su (4.9) 

where 2 denotes summation over all eight nearest neighbors of P and: 
w 


= C a P - C a P + C a P - Cect 1 ? + d B P + d 6 
e e w w n n s s u e e w p 

P 

w 

(4.10) 

(P « f , 0 w U(f w r w , w 

= (C a + d 6 ) + (CcCtc - c a - dj - 
ww ww ss > nn n n 

-stf 

(4.11) 

= (-C a E + d B E ) + (C_a E - C a E - d B E - 
e e e p e ' s u s n n n p n 

d s B s^ 

(4.12) 

3 r NW r nw NW . . NW 

a.,,, = Ca -Ca + d B + dB 
NW w w n n w w n n 


(4.13) 

a r NE f HE, . NE 

a.ir* — "C q — C q dn B— 

NE e e n n n 


(4.14) 


and 

vv 

S = / IS dxdy (4.15) 

u yJ x.J 
J s w 

Several points about the above expressions are noteworthy. Firstly the "a" 
coefficients contain contributions from the fluxes from more than one cell 
face: for example, the a w coefficient contains, in the first brackets, 
contributors from F w and, in the second brackets, quantities relating to the 
intersecting faces s and n. Secondly, there are no obvious bounds on the 
coefficients: they may be positive, negative or zero. The significance of 
these observations will be explained below. 

Having derived the general discretized form of the transport equation, the 
requirements that this equation has to meet to mirror the differential 
equation can be established. These requirements are discussed in the following 
sections. 


16 



4. 1.1.1 Conservation 

Equation 4.1 is, as already noted, an exact statement of the conservation of $ 
principle over the cell volume. Solutions of the discrete version (Equation 
4.9) will obey this principle provided that the weighting factors a and a are 
uniquely defined at the cell boundaries, rather than at grid nodes for the 
reason that only the former practice will ensure that the fluxes are 
continuous across the cell boundaries. 

4. 1.1. 2 Boundedness 

It is well known that, in the absence of sources, solutions to Equation 4.1 
will lie within bounds given by: 

min ( g$ B ) < d <max (dg) (4.16) 

where de denotes the boundary distribution. Sufficient conditions for the 
discrete solution to obey this "boundedness" principle are that: 

o the coefficients a p and aj must all have the same sign, here taken 
without loss of generality to be positive, i.e., 

a p , ai 0 (4.17) 

o the central coefficient must be the sum of its neighbors, i.e., 

a p = z a i (4.18) 

These ensure that when Su = 0, tf p is simply the weighted mean of the 
neighboring so no extraneous extreme can be produced. 

The above conditions lead to further constraints on the weighting factors, 
which may be summarized, by reference to those relating to F w , as follows: 

o they must be non- negative, i.e, 

«w> Bw >_ 0 (4.19) 

o the a w must sum to unity, i.e. 


2 a 
W 


i 

w 


1 


o the b w must obey the resolution: 


(4.20) 


2 _ = 2 + e\ 

w w w w • 


(4.21) 


Inspection of the coefficient definitions (Equations 4.10 - 4.14) reveals that 
obedience to the constraints of Equations 4.19 to 4.21 is not in itself 
sufficient to ensure that the conditions of Equations 4.17 and 4.18 will be 
satisfied. Infringement may occur from either the convection contributions 
(since the convective coefficients are unconstrained as to both magnitude and 
sign, apart from the overall mass continuity requirement) or from the 
diffusion contributions. 


17 



4. 1.1. 3 Transportation 


The directional properties exhibited by fluid transport phenomena are 
well-known and are signaled in Equation 4.1 by the fact that, for S = 0, the 
equations become hyperbolic and streamlines become characteristics as the 
Peclet number, P v L/r , (where L is a typical length scale) tends to 
infinity. The implication for the discretization is that there are domains of 
influence and dependence requirements to be observed at each mesh point which, 
simply stated, entail that at large Peclet numbers information should be 
propagated along streamlines in the direction of the flow. Failure to observe 
this requirement can give rise to unstable solutions (nonphysical 
oscillations) because the coefficients can become negative. 

4.1.2 Accuracy 

Accuracy is the primary requirement of any scheme that proposes to reduce 
numerical diffusion. The scheme should be at least as accurate as the hybrid 
scheme in all situations, and more accurate in most situations. 


The evaluation of the relative accuracies of different schemes is not a 
straight forward task. The ultimate test of a numerical scheme should be based 
on its ability to accurately calculate real flows. In the present context 
these flows should be of the type that occur in a gas turbine combustor. 
Unfortunately, such flows are complex, so analytical solutions are not 
available for comparison. Measurements in these flows are subject to large 
experimental inaccuracies. Moreover, these flows are turbulent, which means 
that the solution accuracy is clouded by the turbulence model being used. 

Traditional means of evaluating accuracy of numerical schemes are available. 
For example, Taylor series analysis will give the order of the accuracy of any 
scheme. Since this analysis is only accurate in the limit when mesh spacing 
tends to zero, it is by no means clear that a higher-order scheme will 
necessarily yield results that are more accurate than a lower-order scheme on 
a mesh of finite size. Comparing solutions against model problems also has 
limited generality because model problems may not be able to represent 
adequately all the characteristics of the class of real flows for which the 
numerical scheme is to be developed. 

The performance of a number of different numerical schemes can be evaluated by 
applying them to a single computational cell for different test cases. This 
method allows different schemes to be assessed quickly against each other, but 
has limited value because the conclusions are valid only for the particular 
test cases studied. 

It is clear that any one approach alone is not sufficient to assess the 
accuracy of proposed schemes, and care should be taken in drawing general 
conclusions from any such single comparison. In this study, all of the above 
methods of accuracy assessment have been used in selecting the most accurate 
scheme. 
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4.1.3 Stability 


Stability of a numerical scheme means that the computations using the scheme 
should be able to converge to a solution. The stability of a numerical scheme 
is, however, also related to the type of solution algorithm being used. For 
example, the present algorithm (an Alternating Direction Impl icit-ADI- method 
using a Tridiagonal Matrix Algorithm-TDMA) requires that the coefficient 
matrix formed from the differencing be diagonally dominant. Any alternative 
differencing scheme for reducing numerical diffusion which generates a 
coefficient matrix in which this dominance is lost will be unstable if the 
present solution algorithm is used, regardless of the overall stability of the 
finite differencing used. However, a more robust solution algorithm in place 
of the ADI method may be able to converge to a solution with the same set of 
difference equations. 

4.1.4 Economy 

Any scheme, if it is to be used in production versions of present computer 
codes, should be economical to use when computing realistic flow problems. 
Economy is understood to mean that the computing time required when using the 
scheme should not be prohibitive. A scheme that yields a marginally more 
accurate solution but requires significantly greater computer time would not 
be cost effective. Computer codes requiring large computing times, usually 
increase the time required to complete a job; consequently, regular usage of 
the code is inhibited. 

4.1.5 Variable Storage Requirements 

Although the memory of present-day computers is extremely large, the number of 
variables required to be stored for a computer program that solves the 
equations for a two-phase turbulent, reacting, radiating flow using 
curvilinear coordinates can be very large, especially if the flow is 
three-dimensional and the geometry is complex. The storage capacity of even 
modern computers can be easily exceeded in such cases. Since the memory 
required is the product of the number of variables stored, the number of nodes 
in the mesh, and the number of nodes in the computational molecule, there is a 
trade-off between the variable storage requirements of a new scheme and the 
increase in accuracy that is achieved by using it. 

With a more accurate scheme, a relatively coarse mesh may be used to achieve 
better accuracy compared to the existing scheme with a fine mesh. Therefore, a 
more accurate scheme may reduce the overall memory requirements even though 
the number of nodes in the computational molecule to be stored may have been 
increased. 

4.1.6 Compatibility With Present Codes 

If a scheme is compatible with the present program it requires less time to 
incorporate it into the code. The present study was limited to developing 
improved finite-difference schemes for TEACH-type codes. Therefore, it was 
essential that the proposed scheme be able to converge with the present 
solution algorithm. 
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4.1.7 Complexity 


A simple scheme takes less time to incorporate into a code than a complex 
scheme. However, this is a one-time expense, and a complex scheme which is 
more accurate, economical and stable than a simple scheme is always to be 
preferred. 

4.2 METHOD OF ASSESSMENT 

It is clear from the above discussion that a differencing scheme has to meet 
several criteria before it can be considered for incorporation into the 
present design codes. It is also clear that traditional Taylor series analyses 
are not sufficient to assess the suitability of a given scheme. In this study, 
both model problem calculations and examination of the properties of the 
difference equations have been used to supplement the Taylor series analyses. 

As explained in the previous section, examination of the weighting factors 
shows whether a certain discretized form is conservative or not. Inspection of 
the coefficients reveal whether the discretized equation satisfies the 
boundedness and transportational properties of the differential equation being 
approximated. Examination of the coefficient matrix gives a good indication 
whether it will be amenable to solution by the ADI method presently used in 
the design codes. As is well known, ADI methods work best when the coefficient 
matrix is diagonally dominant. 

Model problem studies, in addition to assessing the accuracy of a given 
scheme, are useful in estimating its stability, cost effectiveness, complexity 
and variable storage requirements. 

4.2.1 Limitations of Method of Assessment 

The method of assessment described above is more general and informative than 
use of the Taylor series analysis method alone. Two limitations of the method 
of assessment should be noted, however. First, schemes more compatible with 
the present solver (the Alternating Direction Implicit method using a 
Tridiagonal Matrix Algorithm) have an advantage over those schemes that are 
less compatible; in principle, these schemes should be judged independently of 
the solver.. Second, the model problems used herein usually demonstrate the 
performance of a scheme for highly idealized flows possessing only a single 
flow complexity; these simple flows are used because analytical solutions are 
available. However, the practical performance of each scheme can only be 
assessed by computing more realistic flow situations. 

In the next section, difference equations for the candidate schemes will be 
derived and their suitability will be assessed by examining the properties of 
the resulting equations. 
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4.3 DESCRIPTION AND PRELIMINARY SCREENING OF CANDIDATE SCHEMES 

The following candidate schemes were assessed for incorporation into present 
design codes of the TEACH type. 

1. Agarwal Differencing Scheme (ADS) 

2. Quadratic Upwind Differencing Scheme (QUDS) 

3. Skew Upwind Differencing Scheme (SUDS) 

4. Cubic Spline Schemes (CSS) 

5. Glass and Rodi Hermitian Scheme (GRHS) 

6. Flux Blending Schemes 

In the following sections the candidate schemes will be described. The 
descriptions will be brief for the most part, since many of the schemes are 
well known and have been presented elsewhere. The less-familiar flux blending 
approach will be developed in more detail. For the sake of brevity and ease of 
understanding, the equations will be presented for the special circumstances 
of flows in which u, v, p and r are uniform and positive, S=0, and the 
computing mesh is square with spacing a. The evaluation of the candidate 
schemes for their boundedness, transportation, and conservation properties 
will also be discussed in this section. Any improved accuracy scheme must 
satisfy these criteria for it to be considered as a candidate for 
incorporation into computational fluid mechanics codes. 

4.3.1 Central Difference Scheme (CDS) 

By way of illustration, and for later reference, the properties of the 
familiar CDS will be examined in the light of the framework described in 
Section 4.1. 


The non-zero weighting coefficients for this are: 
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(4.22) 


and likewise for a e , B e etc. These give 

a P = sa i 

where: 

a w = C w ( l + F _) 9 a S = C s { 7 + 

e w s 


(4.23) 
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a E = C e { "I + t 4~ ) a N = C n ( " 7 + 
e 

a NW = a SW = a ME = a sE = 0 



In the above equations, P g * C y/ d w - 

c w 

These show that the scheme is conservative, but bounded only for Pe-j < 2. 
The cause of the latter limitation is failure to obey the transportation 
requirement (i.e., use of downstream d values in the convective flux 
approximation) . 


For a one-dimensional problem, dp can be obtained by using equation 4.23. 
The result is: 


1 
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Pe 


4.3.2 Upwind Difference Scheme (UDS) 


(4.24) 


In this scheme the b are identical to the above, but the o are given by the 
following expressions: 



and the coefficients have the form: 
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(4.25) 
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This scheme is seen to be both conservative and unconditionally bounded . Until 
recently it and the familiar Hybrid Differencing Scheme (HDS), which is an 
amalgam of CDS and UDS, were the only known schemes to possess these 
properties. The inaccuracy of these schemes at large Peclet numbers is well 
known, and the schemes that will be described in the following sections are 
intended to improve on this aspect. As will be seen, some of these schemes, in 
attempting to improve the accuracy, lose the essential qualities of the UDS 
and HDS schemes, namely conservation and boundedness. 

4.3.3 Agarwal Differencing Scheme (ADS) 

This scheme, which was developed by Agarwal (Ref. 14 ) employs the extended 
nine-point computational molecule shown in Figure 4-2. The scheme is of the 
finite-difference variety and evaluates the derivatives in the following way. 


Convective Terms - These are expressed via Taylor series as: 


eu l|l t = <„u) h - 6 
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in which (a ) D is approximated in an upwind fashion by: 
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(4.26) 


(4.27) 


(4.28) 


Diffusion Terms 
giving: 


These are evaluated in the usual central difference fashion. 



(4.29) 
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Figure 4-2 Extended Nine-Point Computational Molecule 


Approximation of the y derivatives in a like fashion gives, after assembly, 
the following expressions for the non-zero coeffi cents of the discretized 
transport equation, for u, v> o: 


a E = ‘ 3 C p + d 

a W " C p + d 


*ww 


1 c 

■6 C P 


a = i C n + 2d 

p 2 p 


(4.30) 


where C 

P 


(pu) a 


and d = r /A 


The following are comments on the ADS scheme. 

1. This scheme is non-conservative, since Equation 4.27 does not guarantee 
continuity of convective flux between adjacent pairs of mesh points. 

2. It is also not unconditionally bounded, because certain coefficients are 
either negative or can become so in some circumstances. 
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Concerning the latter point , Table 4-1 gives the form of the coefficients for 
one-dimensional transport in terms of the ratio Pe = c/d, and their limiting 
values (after normalization by an) when this parameter tends to infinity. 

Also given is the absolute sum sa^/ap , which is a measure of the 
diagonal dominance of the coefficient matrix. Similar information is provided 
for the other schemes examined. 


TABLE 4-1 

ONE-DIMENSIONAL COEFFICIENTS AND THEIR LIMITING VALUES FOR INFINITE PE 


LIMIT PE-oo 


SCHEME 

a WW 

a W 

a E 

a EE 

a p 

a WW /a P 

ayy/ap 

ag/ap 

a EE /a P 

E|a|/ap| 

ADS 

Pe 

6 

Pe + 1 


0 

^ + 2 

1 

3 

2 

1 

6 

0 

2 5 

QUDS 

Pe 
“ 8 



0 

3fe + 2 

1 

3 

7 

3 

-1 

0 

3 I 

UDS / 
SUDS 

0 

l(|Pe| + Pe) + 1 

l(|Pe|-Pe)+1 

0 

2 + |Pe| 

0 

1 

0 

0 

1 

CDS/ 

CSS 

0 

5 (<H 


0 

1 

0 

OO 

— 00 

0 

oo 


In the present instance it is clear that the negative coefficients are not 
insignificant in comparison to the positive ones and, moreover, the matrix is 
far from diagonally dominant. The expectation is therefore that problems of 
"nonphysical oscillations" and numerical instability may arise in generating 
solutions for this scheme. Since this scheme is not conservative and 
unbounded, it does not seem to be a promising scheme and was not considered 
further. 

4.3.4 Quadratic Upwind Differencing Scheme (QUDS) 

This scheme, which was proposed by Leonard (Ref. 15 ), is derived via the 
finite volume approach. It approximates the diffusion operator using 
conventional central differencing but uses upstream - biased quadratic 
interpolation to approximate the convective operator. For example, the term 
(pu) w A0 w is evaluated as follows: 


(pu> w A K = (pu) w A 


j’l*WW + ! V'!*p 

( “ 8 **E + 8 + 1 ^p 



(4.31) 
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The presence of the negative coefficients in the above is noteworthy, for it 
illustrates the point made earlier about the possibility of these arising when 
higher-order interpolation is employed. The adoption of these practices for 
the remaining terms gives, for (u,v)ao, the following non-zero coefficients: 


where: 
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i c x + d 


= £ Cy + d 


a WW = ' V 8 

a ss = - Cy/8 
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C x = (pUA) w = ( pUA) e 
C y = ( pUA ) n = (pUA) s 

The following are comments on the QUDS scheme. 


(4.32) 


1. As already noted, this scheme is conservative. 

2. However, it is not unconditionally bounded, as is indicated by the above 
and Table 4-1 which also shows that it departs further from diagonal 
dominance than the ADS in the large Pe limit. Although this scheme is 
unbounded in its present form, it can be bounded using the bounding 
schemes to be described later. Hence this scheme will not be eliminated at 

- present. 


4.3.5 Skew Upwind Differencing Scheme (SUDS) 

This is a finite volume scheme which was developed by Raithby (Ref. 16 ) and 
employs the compact nine-point computational molecule illustrated in Figure 
4-1. In this scheme, the convective flux is obtained by employing upwind 
differencing along streamlines. The streamline direction is defined at each 
boundary by the velocity direction (Figure 4-3). The convected <( is obtained 
by back-projection of this vector until it intersects a grid line and then 
interpolation at the intersection point, the interpolation practice varying 
according to the sector in which the vector lies. By way of illustration, in 
the example of Figure 4-3, if the intersection lies between W and SW, linear 
interpolation is used. Thus: 


< = min (1 , Cy/2 C x ) 
W , SW 

a, = i - a, , 
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Figure 4-3 Control Volume for Skew Upwind Differencing for West Face of 
Control Volume and Positive Velocity Components 


The above practices give rise, for the particular flow directions shown in 
Figure 4-3, to the following non-zero coefficient expressions: 



(4-33) 
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The following comments concern the SUDS scheme. 

1. This scheme is conservative for reasons already explained. 

2. It is not, however, unconditionally bounded because the net convective 
contributions to a w and a s can be of either sign. This condition is a 
consequence of the fluxes leaving through the n and e cell faces being 
linked to and 0 S , respectively. 

3. The scheme reduces to the UDS when the mesh and flow are aligned, so in 
these circumstances it j_s bounded. 

On the basis of the second comment it can be expected that the solutions 
produced by SUDS may exhibit "nonphysical oscillations" and may be difficult 
to obtain due to the matrix not being diagonally dominant. Although this 
scheme is not unconditionally bounded it will be retained for further accuracy 
evaluation because it can be bounded by the bounding schemes to be described 
later. 


4.3.6 Cubic Spline Scheme (CSS) 

The adjective "spline" will be used here to denote schemes based on assumed 
interpolation formulae whose coefficients are linked to nodal values of $ and 
one or more of its derivatives, (s5‘ = m) and («5" = M). A particular scheme 
will be examined briefly, namely the Cubic Spline Scheme (CSS), which has been 
employed by Rubin and Graves (Ref. 17 ), Vacca, Werle and Polak (Ref. 18 ) and 
Kumar (Ref. 19 ), among others. A cubic spline: 


3 

S = 2 a • (x - x > 
i=0 s 


in the interval xj_i < x < Xj, is a polynomial function which is 
continuous in i, m and M sucn that: 


S' 


= M J- 


(x j " 


x) 


(x - 


M j 


x j-l } 


and hence, by integration: 
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(X - Xjj)' 
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x.-x 


-r M j-i t ] ( t ] + u j _M j 3 H 


A 2 x-x . . 

A w J ~ ^ 


(4-34) 


(4-35) 


where a = Xj-Xj_i. The Mj in the above equation are unknown, but can be 
related to the J $j by requiring continuity of Mj at the matching of two 
splines, to give: 
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(4-36) 


M j-1 * 4M J * M JH ■ *2 «J-1 - 2 *j + W 

A 

A similar derivation provides the following link between the mj and 0j, 


m j-i + 4 v vi = a ( Vi - Vi } 


(4-37) 


The above, along with the functional relation between i j , Mj and mj 
obtained by applying the differential equation to be solved at the nodal 
points, i.e.. 


f (0j, m y Mj) = 0 


(4-38) 


produces three coupled sets of equations in the three unknowns. The resulting 
block 3x3 tri diagonal matrix requires in general an inversion procedure for 
solution, so the properties of the equations and the results they will produce 
are not so readily analyzed as were the earlier schemes considered here. 
However, the example application described below is a useful exception to this 
rule. 


We consider the solution by the CSS of the 1-D convection-diffusion equation 
(designated in later sections as test Case 0D1) for which the governing 
differential equation is, in the present notation: 


Pe A 

— - m. - M. = 0 (4-39) 

a J J 

where Pe A = pUA/r is a Peclet number. This equation can be used to eliminate 
the Mj from Equation 4.36, thus: 


eA 


(m j-i * 


4m. 

J 


Vi> 


6 

= IZ 




(4.40) 


The mj can now be eliminated by combining Equations 4.37 and 4.40 to yield: 

Pe A (Vl ' = 0j+i + ^ _ 2 (4-41) 


or, in the "compass" notation and after rearrangement: 

(4 - 42) 

The above equation is identical to Equation 4.24 which was produced by three- 
point centered differencing. 
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The properties of the CDS are well-known and have been reviewed in Section 
4,3,1. Briefly, the scheme is second order, positive and bounded provided that 
Pe A <2. Above this value, diagonal dominance is lost, "nonphysical 
oscillations" occur, and the accuracy generally falls below that of lower 
order schemes. The available evidence about the CSS falls within this pattern; 
it should be noted that, in nearly all applications reported in the literature 
where the scheme has been successfully applied, the cell Peclet number Pe A 
has been at or below the critical level in regions with steep gradients (Ref 
4.5). 


It would appear that the CSS is unlikely to be useful in the 
convection-dominated circumstances of interest in the present study; in some 
cases, it may be disadvantageous. Similar comments are likely to apply to 
higher order centered splines, the common deficiency being infringement of the 
transportation property. For convection-dominated flows, it may be useful to 
introduce the notion of an upwind-biased spline, as has been done in a 
preliminary way by Kumar (Ref. 19), but this is a matter for future research. 
No further consideration will therefore be given to the CSS in the remainder 
of this report. 

4.3.7 Glass and Rodi Hermitian Scheme (GRHS) 


This scheme, developed by Glass and Rodi (Ref. 20), uses a cubic 
interpolation polynomial instead of a cubic spline. This polynomial can be 
written in the following form in two dimensions. 


4 

a> = z 

i = 1 


( Vl + a x,A,i + a y,1 m y,i + \y,i Vi 1 


(4.43) 


where the a's are cubic polynomials in x and y, and m x , my, and M X y 
denote 3 $/ 3 x, 3^/sy and 3^d/3X3y, respectively . This interpolation formula 
is employed in an explicit, time-marching method of characteristics 
calculation of convection whose essential features are listed below. 

1. At each new time level, the characteristic line (i.e. particle path) 
passing through each grid node is traced back to the previous level from 
knowledge of the velocity field. 

2. The value of i at the intersection point of the characteristic with the 
old time level plane is determined using Equation 4.43 for each cell. This 
value is used to evaluate the convective flux. 

3. The values of m* and Mj appearing in Equation 4,43 are 
determined at each mesh point by solving transport equations for them; 
those for Mj being derived through differentiation of the i transport 
equation. 
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The novel character of this method complicates evaluation and comparison with 
the other schemes considered here. Also* this time-marching scheme is not 
compatible with the present steady-state formulation. It is acknowledged by 
its developers to be non-conservative (although they show that if the mesh is 
sufficiently fine to give good accuracy in the conservation errors are also 
small). It is also not unconditionally bounded, but there is no obvious reason 
why the coefficients in Equation 4.43 cannot become negative. In addition, the 
usual stability limits on an explicit method apply, e.g., in one dimension: 



r&t l 

2 < 1 
A 


(4.44) 


The accuracy of the method, for the test cases examined by its developers, is 
markedly superior to those of CDS and UDS with which they were compared. It 
required about twice as much computing time as CDS and UDS on a given mesh, 
but on a time-for-given-accuracy basis it would undoubtedly prove better. 
These, however, were explicit calculations which invariably require 
substantially more time to reach the steady state than do implicit methods. 
The utility of the approach for the present purposes then depends on whether 
it can be implemented implicitly. This matter is beyond the scope of the 
present study, so the GRHS was also excluded from further consideration here. 


4,3.8 Schemes Generated by Flux Blending 

All of the otherwise desirable schemes that have been studied in this section 
suffer from the defect that the coefficients generated by these schemes are 
not unconditionally positive. Negative coefficients are not desirable from 
stability and boundedness considerations. 

Hence, a method is needed that would either make the coefficients uncondition- 
ally positive or at least modify them so that the solution obtained remains 
within bounds. It is apparent that in the absence of this method all of the 
schemes studied so far will have serious flaws. There are only a few such 
methods available in the literature. A. D. Gosman (consultant to this program) 
and his co-workers at Imperial College, London, have been working on what are 
called Flux Blending Schemes. These schemes will be described in detail below. 


4. 3.8.1 The Flux Blending Strategy 

The concept of flux blending was prompted by the "Flux Corrected Transport". 
(FCT) technique developed by Boris and Book (Ref. 21 ) and others, in the 
context of explicit time-marchinq techniques. The approach adopted here for 
implicit steady-state calculations was developed by Gosman and Lai (Ref. 13 ) 
and Gosman and Peric (Ref. 22 ) and differs from FCT in several important 
respects. It operates as follows: in the assembly of a flux-blended scheme, 
the flux F at each cell boundary is expressed as the weighted mean of a 
bounded (but perforce "low order") scheme, F L and an unbounded, but higher 
order scheme, F H , thus: 

F = y F H + (1-y) F L (4.45) 

Here y s the blending factor (0<y<1) is uniquely defined at the cell boundary 
in question, so that the blended scheme is conservative. The other important 
requirement on y is that it should be as large as possible, within the over- 
riding constraint that the solution should be properly bounded. 
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Two alternative strategies have been developed for ensuring that the last 
named requirement is fulfilled (or nearly so). These strategies are listed 
below. 


a) Strategy 1 - This strategy operates by choosing y so as to suppress 
negative coefficients in the blended scheme: this, as was demonstrated 
earlier, is a sufficient condition for a bounded solution. This is the 
approach used by Gosman and Lai (Ref. 13 ). 

b) Strategy 2 - This strategy, which was developed under the present 
contract, ensures that when negative coefficients occur, their 
contribution to the solution is below the level which would cause it to be 
out of bounds. This approach is more akin to the original FCT. 

Each of these strategies has its advantages and disadvantages which are best 
brought out by examining the properties of some specific flux-blended schemes. 
Three of these will now be described, namely BSUDS1, which is a blend of the 
UDS and SUDS schemes based on strategy 1; BSUDS2, which blends the same 
schemes using strategy 2; and BQUDS, which uses the latter strategy to blend 
the UDS and QUDS. 


4. 3. 8. 2 The BSUDS1 Scheme 


It is not difficult to show that if the UDS and SUDS fluxes are blended in the 
manner just described, the resulting non-zero coefficients for (u, v>o) are: 


a,, = d + C (l + y - y ) - y C 
W x ’w w V T n y n 


a £ = d 

a S = d + C y ^ ' T s “s ^ Y e C x “e 
a N = d 

a-,. = y C a SW + y C ct SW 
SW w x w s y s 

a p = 2 a. 


(4.46) 


where the a are defined as before. Note that for the present, 
y e = yw = yx and Yn = Ys = Yy. The above coefficients bear a 
close similarity to those of the parent SUDS scheme, (c.f. Equation 4.33) but 
it is clear that in the present case it is always possible to avoid negative 
coefficients through appropriate specification of the y's. 

The determination of the appropriate y field is however a non-trivial matter, 
for inspection of the expressions for a w and a s (the only coefficients 
which could become negative in the present circumstances) reveals that the 
values are interdependent; i.e. a w and a s contain both y x and y y as do 
certain of the other coefficients. 
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What has in fact emerged is an optimization problem of the following kind. 
Maximize the object function: 


F ( y ) ~ + ^2 

subject to the constraints that: 
0 < y < 1 


* <e i - f i 


.)> 0 


(4.47) 

(4.48) 

(4.49) 


where yi i and y 2 i are the weighting factors appearing in a particular 
coefficient a-j ancl (e-j - f-j yi i) s and g-j are positive quantities 
deducible from the definition of that coefficient. For example, 

e w = d + c x , fy/ = C x ^ ~ > 9w = c y a !f0 

Lai Ref. 23 used the 'Simplex' optimization method for this purpose and 
derived the following expressions for the optimum y. 


r y = rft - 1 

C x 

1, if Pe > 2C /(C - C ) and 1 < ^ < 2 

. T opt . y (4.51) 

x 5 2(1/Pe + 1 - k ) otherwise 

y y 


where: K = min (1, c /2 C ) 

The follow are comments on the BSUDS1 scheme. 

1. This scheme is conservative, has positive coefficients, and is 
unconditionally bounded. 

2. In contrast to the FCT method, it operates only on the coefficients and 
makes no explicit reference to, or requires knowledge of, the solution 
bounds. 

3. Because the resulting coefficient matrix is unconditionally diagonally 
dominant, numerical stability problems are minimized. 

4. It has, however, two disadvantages: firstly, the positive coefficient 
requirement can be argued to be overly- stringent (i.e, it is a 
"sufficient" rather than a "necessary" condition for boundedness) and may 
lead to excessive blending of the low-order scheme; and secondly, it is 
only applicable when the computational molecules of the low and high order 
schemes coincide at those nodes where the latter produces negative 
coefficients. (In this respect the UDS and SUDS schemes are well-matched, 
for SUDS has positive coefficients at all but the principal nodes N, S, E 
and W). These disadvantages are overcome (at the expense of others) by the 
second blending strategy, which will now be described. 
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4. 3,8.3 The BSUDS2 and BQUDS2 Schemes 

As mentioned earlier, Blending Strategy 2 relaxes the positive coefficient 
requirement and replaces it with one of adherence to estimated bounds on <f>. 
This gives rise to two distinct problems; one being to determine the bounds in 
a sufficiently general way, and the other to work out the appropriate y 
values. This latter aspect will be dealt with first, on the assumption that 
the bounds are known. 


The y calculation is iterative and starts from an initial prescription of 
unity values everywhere. In the first stage of each iteration the y values at 
all boundaries of a particular cell are taken to be equal when it is first 
considered. This allows the coefficients of the difference equation to be 
expressed in the following general way in the first iteration: 


a. = a.' + ya." (4.52) 

( a i ' + y a - ' ' ) «J p = a. + y a.. "I. (4.53) 

where a-j 1 contains those terms not depending on y, and ya-{ " contains the 
remainder. Let the solution to Equation 4.53 for given d-j be d p and let 
the known lower and upper bounds, one of which will be imposed^ if d p lies 
outside them, be d mi - n and s5 max respectively. It is then easy to show that 
if: 


a) 


b) 


6 > 6 
p max 


then 


. <l max£ a i' - 

I a . ' ' d . - i y a. ' 
^ i maxi i 


4 < then 
p min 


_ - Wi ■» 

T 


(4.54) 


(4.55) 


if is within bounds, y is simply left unchanged at the prevailing level. 

It should be noted that these equations are valid only in the absence of 
sources and will be modified when sources are present (see Section 5.1.5). 

The foregoing procedure produces a field of y's which is non-unique, for in 
general, two values will have been calculated for each cell face. This is 
resolved in a second stage simply by taking the minimum of the two values in 
the knowledge that this is a "safe" practice, i.e., it will always maintain 
the bounds. This completes the iteration, after which (a 1 -”) m is replaced 
by ( ya-j 1 1 ) m_1 (where m is the iteration counter) in preparation for the 
next cycle. 
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There remains the non-trivial matter of how ^-fn and i^ ax are to be 
determined. This is the subject of continuing research, but quite encouraging 
results have been obtained by requiring that $ should be separately bounded by 
the prevailing neighboring along each grid line passing through P. In the 

model problem studies, this produces two estimates for the bounds, viz., 

( W raax'l - mf "> max ^E* «V (4.56) 

'“min, max’z = l " 1n - ma * ^S 1 (4.57) 

and the practice is to calculate izi via Equations 4.56 and 4.57 for both and 
then take the minimum value. There are other alternatives to this scheme and 
this formulation is the one which was found useful for model problem studies. 
For the two and three dimensional codes, this formulation is replaced by 
another method of computing bounds on 0p (see equations 5-54 and 5-55). 

The application of this blending practice to SUDS and QUDS (for which the UDS 
has been used as the lower order scheme) is now a straightforward matter of 
deducing the appropriate expressions for the a-j ' and a-} 1 '. This has been 
done and is summarized in Table 4-1 I. 


TABLE 4-1 1 

COEFFICIENTS OF BQUDS AND BSUDS2 FOR THE CASE OF 
POSITIVE RADIAL AND AXIAL VELOCITIES 


SCHEME 

COEFFICIENTS (u >0, v>0) 


a E 

a W 

a N 

a s 

a EE 

a WW 

a NN 

a ss 

BQuUo 

d e - g Ye^e 

d w + C w 

” § >w Cw + g 7e c e 

d n “ g 

d s + C s 
~ 8 y s C s + 

0 

- g Yw^W 

0 

1 r 
- g Ys c s 


SCHEME 

COEFFICIENTS (u >0, v>0) 


a E 

a w 

a N 

a S 

a NE 

a NW 

a SE- 

a sw 

BSUDS2 

d e 

d w + C w 

—7wCw ^ — a w ) 
r W 

~ 7 ” C n« n 

d n 

d s + c s 
-Y S c s n-«|) 

- Ye ( - : e a e 

0 

0 

0 

7 wC w aSW + 
+ Y s C s a s S W 
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Comments 


1„ The advantages of this blending procedure are that it is capable in 
principle of producing conservative and bounded solutions for any 
high-order scheme which is itself conservative, i.e. the "overlapping 
molecule" structures of the strategy 1 do not apply. 

2. There are however several disadvantages as compared with strategy 1: 

a) The difference equations are now non-linear (i.e. the a-j's depend 
on the d's through y) even when the parent differential equation is 
not. Therefore, the cost and complexity of obtaining solutions 
increases. On the other hand, in many practical calculations the aj 
are not constants anyway. 

b) The solution is only guaranteed to be bounded when the final result 
is obtained, so precautions must be taken to avoid the adverse 
consequences of overshoots at intermediate stages (see Appendix A-5). 

c) The coefficient matrix will not necessarily be diagonally dominant 
(although whatever lower-order blending is performed will increase 
diagonal dominance) so numerical stability problems may still arise. 

d) The uniqueness of the solution cannot be guaranteed due to the 
non-linearity of the equations and somewhat arbitrary although 
logical method of estimating the bounds. Studies thus far suggest 
that this is not a serious problem in practice. 


Thus far, attention has been focused on the conservation and boundedness 
properties of the schemes chosen for evaluation. Based on this preliminary 
screening, three schemes, namely the ADS, the CSS and GRHS, were excluded from 
further evaluation; possible defects (negative coefficients) were identified 
in certain of the others. Recognition of these defects led to the development 
of flux blending schemes and these have been described in detail. 

4.4 ACCURACY EVALUATION 

This section contains a discussion of the accuracy evaluation of the remaining 
schemes. The consequences of the defects identified in the schemes in the 
previous section will also become apparent. The flux blending strategies will 
also be applied to some of the chosen schemes and their effect on the 
accuracy, stability and boundedness of these schemes will be studied by 
calculating some test cases. 

4.4.1 Truncation Error Analysis 

The formal Taylor Series Truncation Errors (TSTEs) of the schemes chosen in 
the previous section for further evaluation will now be presented and 
discussed. Since it is the accuracy of representation of the convective terms 
which is the critical issue, these TSTEs have been evaluated for the limiting 
case of purely convective transport. To facilitate interpretation of the 
results in terms of their numerical diffusion implications, the flow has been 
chosen to be uniform and parallel. 
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The truncation error expressions are given in different forms in Tables 4-1 1 1 
and 4-IV. For obvious reasons no entries appear for schemes based on the flux 
blending strategy 2, i.e. BSUDS2 and BQUDS2 s but it can be assumed that the 
associated TSTEs will be strongly biased towards the higher-order component of 
the blend. Evidence to this effect can be seen in the TSTE for BSUDS1 , 
discussed below. Although ADS has been eliminated from consideration, it is 
included in this discussion because it turns out to be the most accurate 
scheme based on TSTE analyses and reinforces the argument against using only 
Taylor series analyses for accuracy evaluation. 


TABLE 4-1 1 1 

TRUNCATION ERRORS FOR THE CONVECTIVE TERMS OF THE UDS, SUDS, BSUDS1 , QUDS AND 
ADS FOR A UNIFORM FLOW AND SQUARE MESH: CARTESIAN COORDINATE VERSIONS 


SCHEME 

TRUNCATION ERROR 

UDS 

-pA j 
2 

^ d 2 <t> 3 2 <A 

u— f + v— § 
y dx 2 3y 2 J 

SUDS 

-pA 

2 

d 2 <t> d 2 <t> „ 

u— £ + v— ^ + 2v 

3x 2 dy 2 

( \ , , \ d 2 <t> , . /, u \ 

- + k) ;k = min(1, — I 

\2 / 3x3yJ \ 2v J 

BSUDSl 

IpA 

d 2 (/> d 2 <t> „ d 2 <t> ~| 

dx^ dy z dxdy J 

QUDS 

~ — pA^ 

24 

3 3 </> d 3 <t> 

u — = + v — = 
3x3 3y3_ 


ADS 

-1 P a3 

12 

3 4 <*> , 3 4 <*>~ 

3x^ 3y4 



The entries of Table 4- 1 1 1 express the TSTE's in the original Cartesian 
coordinates (x, y) and the corresponding velocity components (u, v) in which 
they were derived. However it is more instructive to transform them into the 
"streamline" coordinates (s,n), where s and n are measured along and normal to 
the streamlines, respectively. Table 4-IV contains the transformed versions, 
certain terms of which are singled out for attention by enclosing them in 
square brackets. These forms involve derivatives in the streamline direction 
which, in the present convection-dominated circumstances, will be negligible 
in the absence of sources, thus, it can be argued, it is the remaining terms, 
involving cross-stream derivatives, which are the principal source of error, 
(provided of course that the TSTE is truly representative of the error, a 
point which will be discussed later). 
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According to Table 4- IV the hierarchy of the various schemes is as follows: 1) 
UDS, SUDS and BSUDSl-lst order; 2) QUDS - second order; 3) ADS- third order, 
where the "order" refers to the exponent on the mesh spacing a. If, as in many 
finite difference analyses "order" is erroneously taken to be synonymous with 
accuracy", then the ADS would stand out as superior to the other schemes and 
the UDS, SUDS and BSUDS1 would be collectively dismissed as having the same 
low accuracy. This interpretation is now widely recognized as being misleading 
and evidence to support this view will be given here. 


TABLE 4- IV 

TRUNCATION EXPRESSIONS OF TABLE 4-III, TRANSFORMED INTO (s, n) 

STREAMLINE COORDINATES 


SCHEME 

TRUNCATION ERROR 

UDS 

1 i \/2 02 A T 02 A 02 a “]) 

P|v|A< sln(20) sin(n/4+0) + (cos 3 0+sin 3 0) 1 - sin(20) (sin0 -cos0) 

2 V 2 3 n 2 L g s 2 3s3n 

SUDS 

-i- P |V|A | _1_ sin (20)(cos0— 2sin(0K) ■£_!. rcos 3 0+sin 3 0+ sin0sin(20) ( 1 +K) ^ ^ 

2(2 3n 2 L 3s2 

+ (sin(20) (sin0-co*0H-2sin0 cos(20)( ^+K) AJL j > : K = min (1, 1 -H-l 

3 s 3 n J ) 

BSUDS1 

1 (V 2 " 3^$ 

■ — PlvlA < sin(20) sin (7T/4 - 0 ) 

2 ' ( 2 3n 2 

/ \ 3 2 <l> 

f (sln( 20 ) (sin 0 cos 0 )+ 2 sin 0 cos( 20 ) j ^ 

r 3 2 $ 

+ (cos 3 0 +sin 3 0 +sin 0 s!n( 20 ) 

^ L 3s2 

QUDS 

PlvlA 9 I 4 - -«ni 1 . 404 . 4oi 2 • _ i.fli | ® ,;r.2/9fl) ^ t. 

94 P l v | A* sin{40) ~~L + Icos^e+sin W sin (4B) + sm^oi 

** j 4 • a n 3 [_ * 3s 3 2 3s 2 3n 2, d sSn^J 

ADS 

( g4^> 

— P 1 V / A 3 < -L sin 20 (sin 3 0 + cos 3 0) ^4 + 

+ | sin 2 ( 20 )cos 0 + sln 0 ) g s 2 g n 2 + sln2 12 

(cos 3 0 + sin 3 0 )'T - J + 2 sin( 20 )(sin 3 0 — cos 3 0 ) r — 

a *V 3s 3 3n 

3 4 ^ 1 

0 ) (cos 0 — sm 0 ) > 

3s 3n 3 J) 


It is conventional to interpret those TSTE's which contain second derivatives, 
e.g. a 2 ?l/ 3 x 2 , a 2 ?i/ay 2 , a 2 tf/an 2 as being diffusion-like, by analogy 
with the "real" diffusion terms in the full transport Equation 4.1. This 
analogy allows a numerical diffusion coefficient r num to be defined as the 
rati o : 


r 


num 


= TSTE/(a 2 d/ax 2 ) 


(4.59) 














and evaluated for schemes having TSTE's of this kind. These comprise the 
first-order UDS, SUDS and BSUDS1. Table 4-V gives the numerical diffusivity 
formulae, which are all of the general form 

r num = p A ^ f(e) (4.59) 


TABLE 4-V 

NUMERICAL DIFFUSIVITY EXPRESSIONS FOR UDS, SUDS AND 
BSUDS1 IN STREAMLINE COORDINATES 


SCHEME 

r 

1 num 

UDS 

sj2 — 

-^-pA |V|sin(20)sin(x/4 + 6) 

SUDS 

^-pA |V|sin(20)[cos(0) - 2sin(0k)I; k = min ^1 

BSUDS1 

pA |V|sin(20)sin(7r/4 - 0) 


where |v| = (u 2 + v 2 )!/2 and e = tan -1 (u/v) are, respectively, the 
magnitude of the resultant velocity and the angle which it makes with the 
mesh. The function f(e) is dimensionless and is plotted in o Figure 4-4 for each 
of the schemes mentioned. The plot is terminated at d = 45° because the 
functions are all symmetrical about this angle. It can be seen that i) the 
levels and variations ofr nU m produced by the three schemes are quite 
different, indicating that their "order" is an incomplete indicator of their 
relative performance; ii) the SUDS is substantially better than the UDS, and 
BSUDS1 is between the two, but strongly biased towards the former; and iii) 
whereas rnum f° r t* 16 UDS increases monotonically from zero at e = 0 to a 
maximum at © = 45°, the other schemes have their maxima at intermediate angles 
and have zeros at the two extremes. On the basis of these results and the mode 
of operation of BSUDS2 it can be conjectured that it would probably fall 
between SUDS and BSUDS1. 
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Figure 4-4 Comparison of Numerical Diffusivities for Discretization Schemes 


As a final comment, it should be noted that, while the TSTE's of the ADS and 
QUDS are not diffusion-like and therefore do not permit the evaluation of a 
r num» it is not correct to imply that they will not produce numerical 
smearing or other unfavorable effects; this will be shown in the test cases 
described in Section 4.4.3. 

4.4.2 Description of Test Cases for Accuracy Assessment 

The test cases chosen here to evaluate the discretization schemes are ones for 
which it is possible to obtain exact solutions to the governing equations. The 
cases fall into two classes according to whether they are one- or 
two-dimensional. Within the latter there is a further sub- division according 
to whether "single-cell" or "field" calculations are performed. In the 
first-named kind, a single computational molecule is considered and is 
analyzed using the difference equation by inserting the exact analytical 
values for the neighboring d and comparing the result with the exact i. The 
fact that an algebraic formula is obtained for the $ error is convenient for 
analysis, but the approach has the disadvantage of not showing the error 
propagation and accumulation effects that often occur in practical 
calculations. Hence, "field" tests, in which calculations are performed over a 
representative computing mesh with the exact values imposed at the boundaries 
only, were also included. 


40 



In this section, all of the test cases are described. Results of the model 
problem studies based on these test cases are discussed in Section 4,4.3. 
Unless otherwise stated, these test cases pertain to situations in which P) u, 
v and are everywhere uniform (which means that the flow is everywhere 
uniform and parallel). The computing mesh has square cells of side a. It 
should be noted that most schemes were not applied to all test cases: Table 
4- VI shows the various combinations examined. 


TABLE 4-VI 

SUMMARY OF THE COMBINATIONS OF DISCRETIZATION SCHEME AND 

TEST CASE EXAMINED 


Model problem 


Scheme 

OD1 

OD2 

OD3 

TDF1 

TDF2 

TDF3 

TDF4 

TDS1 

TDS2 

QUDS 



n 



iS 

iS 



SUDS 





V* 

iS 

iS 

V* 

iS 

BSUDS1 




iS 


is 


iS 


BSUDS2 




V* 






BQUDS 


tS 


tS 






UDS 




iS 

y* 




iS 


4. 4. 2.1 One-Dimensional Cases 

These represent solutions to the one-dimensional transport equation. 


drf 1 d 2 $J 

— * " Ve 

dx e L dx* 



(4.60) 


where x* = x/L, Pe L = P uL/r , S* = SLpu and L is a characteristic length. 
The following test cases were examined: 

(a) Convection and diffusion in absence of source. Test Case (0D1). 
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Here S* = 0 and the boundary conditions d = 0, 1 at x* = 0, 1 are imposed to 
give: 


exp (Pe L x*) - 1 
^ “ exp (Pe L ) -1 


(4.61) 


(b) Convection with no diffusion and linear source, Test Case (0D2). 

This is similar to case 0D1 above, but the source S* is now a pi ecewise-1 inear 
distribution. This test case was also employed by Leonard (Ref. 15 ) in his 
presentation of the QUDS scheme. The exact solution is of the form; 


( = d Q + a (x - x Q ) + b (x - x Q ) 2 , for x > x Q (4.62) 

where x* is the starting point of the source and a and b are constants linked 
to its level and gradient. 

(c) Convection with no diffusion and plane source (0D3). 

For this case Pe is assumed to be infinitely large, S*=l at x*=0.2 and is zero 
elsewhere, and d = 0 at x*=0. The solution is: 

d = 0 , x* < 0.2 
d = 1 , x* > 0.2 

4.4.2. 2 Two-Dimensional Cases 

For these test cases, solutions are obtained to the two-dimensional transport 
equation. 


pu 
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3X 


+ -^f) 

ay 


(4.63) 


for the different circumstances described below. These cases will be denoted 
by "TD" with a suffix S or F to indicate whether they were used for single 
cell or field evaluations, respectively. 

(a) Convection in Absence of Diffusion 

Three cases are considered for the convection of a prescribed distribution d 
(n), where n is the coordinate normal to the streamlines. In two of the cases 
the flow is uniform and parallel, Figure 4-5, but in the third case, the 
streamline pattern corresponds to irrotational plane stagnation flow (Figure 
4-6) in which the d are found to be rectangular hyperbolae, i.e. 

d = xy (4.64) 


42 





In all cases the exact solution is simply 3«S/aS = 0 where S is the streamline 
coordinate: thus on each streamline d is constant and equal to its initial 
val ue. 

(i) Single Cell Calculation of Convection of a Gaussian Distribution in 
Uniform Flow (TDS1) - Here j5(n) = exp (-1/2 (nA*) z ) where a* = a /g, 
a is the standard deviation of the Gaussian, and n is measured from 
the streamline passing through node P, Figure 4-7, 

(ii) Field Calculation of Convection of a Step Function in Uniform Flow 
(TDF1) - In this case 


d(n) =0 , for n < n s 
1 , for n > n s 

where n s denotes a selected reference streamline. 


(4.65) 
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(Ill) 


Field Calculation of Convection of a Square Wave Function in 
Irrotational Plane Stagnation Flow (TDF2) 


d(n) = 1 for n s i < n < n S 2 


(4.66) 


= 0 otherwise 

where n s i and n s 2 denote reference streamlines. 

(b) Field Calculations of Convection and Diffusion of Initial Step Profile in 
Uniform Flow (TDF3) 



Figure 4-7 Single Cell Calculation of Convection of a Gaussian Distribution 


This case is similar to TDF1 but now cross-stream diffusion is allowed to 
modify the initial step profile as it proceeds downstream. The exact solution 

is: 

tfU, n) = [1 + erf [-^ n (Pe/a)^]] 

where erf is the error function. 

(c) Scalar Transport in Jet Flow Field (TDF4) 

In this particular case the flow field is non-uniform, being that of a self- 
preserving laminar plane jet. The associated velocity field and the 
corresponding self-preserving scalar field, for the case where the jet is 
initially "heated", are given by Schlicting (Ref. 24) and will not be quoted 
here. In the test calculations the exact velocities are imposed throughout the 
grid in such a way as to ensure that the discrete mass conservation 
requirement is obeyed at each cell; the exact ^values are imposed at the 
boundaries in the usual fashion. 
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(d) Single Cell Calculation of Convection and Diffusion of ^with Initial 
Linear Distributions (TDS2) 

In this case the linear t distributions are imposed at the boundaries. The 
solution for <£ p is given by Castro (Ref. 25) and is in the form of an 
infinite series. 

4.4.3 Results of Test Cases 

The results of some of the test calculations will now be presented under the 
headings of the different schemes which have been tested. In many instances 
the predictions of the UDS are included to provide an indication of the 
improvements which have been achieved. 

4. 4. 3. 2 The QUDS and SUDS 

The QUDS results for the single-cell test Case TDS1 (convection of a Gaussian 
distribution) will be examined first: they are displayed in Figure 4-8a, 
showing 0 p plotted against the angle between the flow and the grid, for 
various values of the ratio a* of the mesh spacing to the standard deviation. 

Evidently at any given a* there is agreement with the exact solution (given by 
gjp = 1 9 for all o e) at e = 0, but a finite and increasing error emerges for 
angles up to 45°, indicating that numerical smearing is occurring. The UDS 
results of Figure 4-8b show similar behavior, but the errors are considerably 
larger. The single-cell test case results for SUDS are shown in Figure 4-8c. 
This figure indicates maximum error levels comparable to those obtained using 
QUDS, but they occur for SUDS at angles of around 15°. This is in contrast to 
QUDS which gives large errors at the larger angles. The latter behavior is in 
conformity with the numerical diffusivity characteristics displayed earlier in 
Figure 4-4. 



(b) UDS 


(c) SUDS 





6 degrees 


0 degrees 


AL 

\ 0.5 

j 0.75 

1 1.0 

Exact solution 


Figure 4-8 Comparison of QUDS, UDS and SUDS for Test Case TDS1 
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Figure 4-9 shows a plot of the error with a* and indicates that the error 
decreases more rapidly for the QUDS than it does for the UDS as the mesh is 
refined. The slopes of the lines are 3 and 1, respectively, in conformity with 
the indications of the TSTE analysis of Table 4- 1 1 1 for the UDS but not 
apparently for QUDS which the analysts indicates to be second order; for a 
Gaussian distribution, the leading error term for QUDS is zero (i.e. 
ad 3 /ax 3 = 0) whereas the next highest term in the error series is non 
zero. Note also that the line representing SUDS in Figure 4-9, has the same 
slope as that of QUDS (the fact that it lies below the QUDS result is of 
lesser significance, since it is merely due to the SUDS error being small at 
the flow angle 30° for which the results are plotted) indicating that in the 
present instance SUDS behaves like a third order scheme. This apparent 
conflict with the truncation error analysis can be explained by the fact that 
leading terms of the error series are in the present case small in comparison 
to the higher order ones. The results of Figure 4-9 show the inadequacy of 
using truncation error analysis alone in determining the accuracy of various 
schemes. 


HDS or UDS 



Figure 4-9 Error as a Function of Mesh Spacing for Test Case TDS1 for 
Various Schemes for Flow Angle of 30° 
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For the field calculations. Figure 4-10 shows the QUDS predictions and exact 
solutions for the convection and diffusion of a step profile (Case TDF3) for 
an overall Peclet number of 50 and for various e. The results are displayed in 
terms of the 0 profiles along the vertical mid-plane of a 9 x 9 grid and show 
the much-reduced smearing error of QUDS as compared with UDS as e is increased 
(it is nonetheless finite and increases with e). Also shown are the 
undesirable undershoots and overshoots of the higher-order QUDS scheme. 
Examination of the SUDS field calculations for Case TDF3, also shown in Figure 
4-10, confirms the general characteristics already described but also 
illustrates the fact that, like QUDS, this scheme can produce unbounded 
behavior when d gradients are steep: in this instance the peak excursions of 
the two schemes are comparable. QUDS would better preserve the slope of the 
exact solution at angles up to 30° but beyond this SUDS is clearly superior. 
Indeed, a general trend can be noted in which a given scheme either produces a 
bounded solution of reduced slope or an unbounded solution of correct slope. 




Figure 4-10 Comparison of QUDS and SUDS Schemes for Test Case TDF3 


A similar evaluation applies to the results for infinite Peclet number (TDF1) 
when SUDS results (Figure 4-11) are compared with the QUDS results (Figures 
4-12). However, at larger angles the SUDS profiles exhibited more severe 
oscillations than those generated by QUDS. 
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The comparisons for the convection of a square wave in a irrotational 
stagnation flow field (TDF2) are shown in Figure 4-13 in terms of the 4 
profile across the outlet boundary of the calculation domain, which is 
identical with the inlet distribution. The QUDS predictions are considerably 
better than those of the UDS but exhibit significant undershoots and 
overshoots. The SUDS predictions for Case TDF2 can also be seen in Figure 4-13 
and compared with the QUDS results in the diagram. On balance there is little 
to recommend either scheme: both reproduce the square wave to a similar 
extent, albeit without maintaining the proper bounds. SUDS maintains the wave 
amplitude better while introducing a phase error (the prediction is shifted to 
the left relative to the exact solution) whereas the reverse is true with the 
QUDS. 



Irrotational stagnation point flow — Scalar profile transport 


Figure 4-13 Comparison of QUDS and SUDS Schemes for Test Case TDF2 


The results obtained using QUDS and SUDS to compute the scalar transport in a 
laminar jet flow (TDF4) are shown in Figure 4-14, where the profiles along the 
vertical mid-plane are plotted. In this case the QUDS solutions remain within 
expected bounds (probably due to the smoothness of the 4 profile) and are more 
accurate than those of the UDS at both angles; however, QUDS exhibits 
substantial diffusion-like errors at the largest angle (e = 45°). SUDS, like 
QUDS, maintains the proper bounds because the 4 gradients are less severe. In 
familiar pattern, QUDS is better at small angles and SUDS is better at large 
angles. 
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Convection and diffusion of a Scalar distribution 


Figure 4-14 Comparison of QUDS and SUDS Schemes for Test Case TDF4 


The comparison of QUDS and SUDS can be summarized in the following manner: 
both schemes suffer from boundedness problems when conditions with regard to 
the cell Peclet number, flow angle relative to the mesh, and gradients in d 
are large. QUDS is more accurate for smaller flow angles while SUDS is more 
accurate for larger flow angles. 

4.4.3. 3 Flux-Blended Schemes 

In this section the performance is examined for BSUDS1 (which was derived by 
blending UDS and SUDS to maintain positive coefficients) and for BSUDS2 and 
BQUDS2, (where the blending criterion is the preservation of local bounds on 
the solution). 

(a) One-Dimensional Tests 

Since in one-dimensional flow, the UDS and SUD are identical, calculations 
were performed only using BQUDS2 and the results are displayed in Figures 4-15 
(Case 0D1) and 4-16 (Case 0D2). In both cases, it can be seen that 
oscillations have been suppressed and the proper bounds maintained; moreover, 
these results were obtained with steep gradients in the dependent variable. 

(b) Two-Dimensional Tests 

A comparison of the BSUDS2 and BQUDS2 solutions for TDF1 of Figure 4-17 shows 
that both exhibit similar characteristics as those exhibited by BQUDS2 for the 
one-dimensional test cases. The tendency in all cases is to maintain the 
proper bounds without introducing excessive numerical diffusion; this is 
indicated by the fact that in flow regimes where the original QUDS and SUDS 
solution were within bounds (Figures 4-11 and 4-12), those solutions are 
virtually unchanged except from some "rounding" of the profile corners. Again 
BQUDS2 seems to be more accurate at small angles whereas BSUDS2 is more 
accurate at large ones. 
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(a) (b) 



• Exact 

— □ — QUDS 

- “®" BQUDS 

~ Denotes solution is 
out of bounds 


Figure 4-15 Comparison of Bounded and Unbounded QUDS Schemes for 
Test Case 0D1 



Pure convection of a Scalar with plane source in one-dimensional flow 


Figure 4-16 Comparison of Bounded and Unbounded QUDS Schemes for 
Test Case 0D2 
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Y 

Pure convection of Scalar profile in a uniform flow field 



Pure convection of a Scalar profile in a uniform flow field 


Figure 4-17 Comparison of BSUDS2 and BQUDS2 Schemes on Test Case TDF1 


The performance of BSUDS2 and BQUDS2 for model problem TDF2, Figure 4-18, is 
similar to their unbounded counterparts. Figure 4-13, except that the 
unrealistic over and undershoots have been suppressed. 



Irrotational stagnation point flow — Scalar profile transport 


Figure 4-18 Comparison of BSUDS2 and BQUDS2 Schemes on Test Case TDF2 


The performance of the alternative blending approach, as embodied in BSUDS1, 
will now be examined. In this instance results are available for the two 
single cell cases and are included for Case TDS1 in Figure 4-9 and the 
accompanying error plot of Figure 4-10. According to the first two figures the 
BSUDS1 errors lie between those of the UDS and SUDS, with a bias in Figure 4-9 
towards the latter as expected from the numerical diffusivity analysis (Figure 
4-1). Results for TDS2 are similar and are not shown. 


The field calculations for Case TDF3 appear in Figure 4-19 and confirm that 
BSUDS1 does not exhibit the SUDS overshoots. It has been found that the 
BSUDS1 solution tends to lie closer to SUDS than to UDS, especially at the 
larger e, due to the fact that the negative coefficients of SLins reach their 
peak magnitudes at the smaller angles. 

The BSUDS1 calculations for the jet case (TDF4) they differed very little from 
those of SUDS at all angles, presumably because the negative coefficients of 
the latter contribute little to the solution due to the smaller gradients 
there. These calculations are not shown. 
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Convection and diffusion of Scalar distribution — initial distribution step 


Figure 4-19 Comparison of SUDS and BSUDS1 Schemes on Test Case TDF3 


It is concluded that the bounding schemes can suppress overshoots and 
undershoots successfully. The loss in accuracy, or increase in numerical 
diffusion, for Bounding Scheme 2 seems to be negligible for the test cases 
considered. Bounding Scheme 1 does introduce numerical diffusion but it is 
still far more accurate than UDS. BQUDS2 seems to be more accurate at small 
flow angles whereas BSUDS2 seems to be more accurate at large flow angles. 

4.5 STABILITY 

Within the context of the present tests, no instability problems were 
encountered with the SUDS and BSUDS1 schemes. For the latter, this behavior is 
expected since the associated coefficient matrix is always diagonally 
dominant. Although the SUDS matrix does not possess this property, the 
departures were evidently not so great as to prevent convergence. However, it 
should be noted that other users have reported such problems. Some problems 
were overcome by using better grids and modifying under-relaxation factors. 

Convergence problems were encountered during the development of Bounding 
Scheme 2. These instabilities were traced to the failure to insure that, for a 
given set of blended coefficients, a reasonable approximation is obtained for 
the corresponding $ field; otherwise, the method of estimating the local 
solution bounds cannot operate properly. This problem, which should really be 
attributed to the properties of the equation solver employed rather than to 
the discretization method, was solved by arranging an inner iteration cycle in 
which, for given coefficients, the residual error of the 6 solution was always 
reduced to a preset level. Once this was done, no further problems were 
encountered. 
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The calculations using the QUDS scheme were prone to instability, whether the 
Tri diagonal Matrix Algorithm (TDMA) or the Pentadi agonal Matrix Algorithm 
(PDMA) were used in the ADI solution procedure. No such problems were 
encountered in ID calculations with PDMA since it then acts as a direct 
solver. It was found that the instabilities could be suppressed through the 
use of conventional under-relaxation (typical factors being 0.5-0. 6) but, as 
will be seen in the next section, the rather large number of iterations 
required suggests that stability is still marginal. BQUDS2 also suffered from 
this problem, but to a slightly lesser extent. 

An important implication of the foregoing findings is that better solution 
algorithms are needed for QUDS-type schemes to make them fully cost-effective. 
It was hoped that the PDMA might provide the answer, since it uses for all of 
the 0 values along a grid line to which it is applied. That this expectation 
was not realized may be due either to the fact that linkages also exist on two 
neighboring lines on each side or to the occurrence of negative coefficients 
at the principal and outlying nodes. It is interesting to note that the more 
stable SUDS generates a more compact 9-diagonal matrix; although negative , 
coefficients may occur, these occur only at the principal nodes. 

4.6 COST EFFECTIVENESS 

The measure of cost effectiveness of the schemes used here is examined in 
terms of their computing time and storage requirements, however, since current 
trends in computer technology indicate that increased storage is becoming more 
plentiful and inexpensive, the emphasis in the cost-effectiveness calculation 
is on computing time. 

Strictly speaking, the comparison between schemes should be on an equal 
accuracy basis, i.e. tests should be performed on each scheme with different 
grids to obtain a pre-specified level of error. However because this is time- 
consuming and also because the differences between the predictions of the 
"higher order" schemes on a given grid are generally much smaller than those 
between any other scheme and UDS, the results for calculations using the same 
grid will be compared. This comparison gives only a first estimate of the cost 
of using alternative schemes regardless of accuracy. Additional cost 
information is provided for the schemes incorporated into the 2D and 3D 
versions of TEACH for the more realistic test cases discussed in Section 6. 

4.6.1 Computing Times 

A representative set of computing times is shown in Table 4-VII. These pertain 
to calculations for Case TDF2 on the 9x9 grid, and were obtained on a CDC 
Cyber 174 machine: however it is not the absolute times which are of interest, 
but rather their relation to the baseline figures for UDS which are also shown. 

Table 4-VII shows the following: (i) the time required for a single assembly 
of all coefficients, T coe ff; (ii) the time required for a single application 
of the ADI solver, T S0 -] V p r , this turns out to be almost the same for both 
TDMA and PDMA versions; (iii) the total number of outer iterations performed, 
^Quter* which is also equal to the total number of coefficient assemblies; 

(iv) the total number of inner iterations performed, N.j nner , which is equal 
to the total number of ADI passes; and finally (v) the total CPU time, 
excluding compilation and initialization. It should be noted that N-fnner and 
N ou ter are identical for all but the BSUDS2 and BQUDS2 schemes; for these 
two, it was necessary to ensure that the inner loop be executed to a preset 
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residual tolerance rather than for a preset number of ADI sweeps (specified as 
unity for all other schemes). 


TABLE 4— V 1 1 

COMPARISON OF COMPUTING TIMES AND NUMBER OF ITERATIONS 
REQUIRED FOR CALCULATIONS OF CASE TDF2 ON 9 X 9 MESH 


Scheme 

Coeff 

assembly 

time 

Solver 

time 

No. of 
outer-iter 

No. of 
inner-iter 

CPU 

time 

UDS 

1.0 

1.0 

1.0 

1.0 

1.0 

SUDS 

7.0 

1.0 

1.0 

1.0 

1.5 

QUDS 

6.0 

1.0 

14 - 20 

14 - 20 

10 - 15 

BSUDS1 

8.0 

1.0 

1 

1 

1.5 

BSUDS2 

8.0 

1.0 

2-3 

5-6 

4 - 5 

BQUDS2 

7.5 

1.0 

7-13 

65 - 160 

32-69 


• Inefficiency of BQUDS2 is due to the solver (TDMA) 
which is not compatable 


Some interesting general features emerge from this table. Firstly, the 
variations between the coefficient assembly times for "higher-order" schemes 
are small in relation to the increase over the UDS figure, the ratio over the 
latter being in the range 6-8. Secondly, T CO p-ff and T so i ver are roughly 
the same, excluding those for the UDS; and thirdly, QUDS and BQUDS2 clearly 
require substantially more outer iterations than the others. 

Turning now to detailed comparisons, it is clear that SUDS and BSUDS1 are the 
most cost-effective for the present applications due to the smaller number of 
iterations (inner and outer) required; the CPU time increase over that for UDS 
is only 50 percent. The additional cost of the first flux-blending procedure 
is evidently insignificant. 
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The next most cost-effective method is BSUDS2, which requires 3 times the 
SUDS/BSUDS1 CPU time, due to the increased number of iterations required for 
convergence; however, BSUDS2 is measurably more accurate than SUDS or BSUDS2. 

QUDS and BQUDS2 are considerably more expensive, requiring respectively 2-3 
and 6-14 times the BSUDS2 CPU time. The additional expenditure is due to the 
increased values of Mouter- f r ° r BQUDS2, there is the additional cost due to 
the substantial level of N-j nner . 

It should be emphasized that, although the results in Table 4- VI I are 
representative of the present tests, the poor performance for QUDS, B0UDS2 and 
BSUDS2 is due in large measure to the unsuitability of the solvers employed. 


TABLE 4- VI 1 1 

SUMMARY OF ADDITIONAL STORAGE REQUIREMENTS 
OVER THOSE OF THE UDS 


SCHEME 

NO. OF ADDITIONAL 
2D ARRAYS 

QUANTITIES STORED 

SUDS 

4 

a NW' a SW' a NE' a SE 

QUDS 

4 

a WW- a SS' a NN' a EE 

BSUDS1&2 

BQUDS2 

6 

AS ABOVE, PLUS 7e/w AND 7n/s 


4.6.2 Computer Storage 

For 2D applications, none of the schemes examined generates unacceptable 
additional storage requirements over those of UDS. Table 4-VIII indicates the 
number of additional two-dimensional arrays required. In the case of SUDS and 
QUDS, the increase is equal to the number of additional "neighbor" 
coefficients, i.e., 4. The flux blended versions require two additional arrays 
to store the 7 values for each coordinate direction. Assuming that there is a 
choice between savings in CPU time, code complexity and savings in storage, 
savings in CPU time, code complexity will take precedence. 
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4.7 COMPLEXITY AND COMPATIBILITY WITH EXISTING TEACH METHODOLOGY 

The evaluation of complexity is necessarily subjective, but if the schemes 
examined were so ranked, they might appear as follows, with complexity 
increasing from left to right: 

ADS 

UDS, SUDS , BSUDS1, BSUDS2 
QUDS BQUDS2 

In no case can it be said that any particular scheme is difficult to 
comprehend and implement, although the possibilities for error are obviously 
greater as one moves away from the simplicity of UDS. Perhaps the major 
pitfalls lie in the imposition of the boundary conditions, which require 
particular care with schemes like QUDS, SUDS and their bounded counter parts. 
This difficulty is common to all higher-order schemes. 

Concerning compatibility with TEACH, it is sufficient to note that all of the 
tests reported herein were performed with versions of this code employing the 
same basic methodology. Some minor problems were encountered, such as those 
due to the extra row of external nodes at boundaries when QUDS and BQUDS2 were 
implemented, but these were more tedious than difficult. Of course, the 
solvers which have been employed are not particularly well suited to QUDS, but 
the modular structure of the code readily allows a better solver to be used. 


4.8 PROGNOSTICATIONS FOR THREE DIMENSIONAL CALCULATIONS AND OTHER CONSIDER- 
ATIONS 

Before the final evaluation is made it is important to consider two further 
aspects which have not been examined in detail in the studies described thus 
far, namely: 

(a) The suitability of the various schemes for three-dimensional applications. 

(b) Experience of application of the schemes to flow field calculations. 

Both the SUDS and QUDS (and bounded versions thereof) have computational 
"molecules" of equal size in 2D, i.e. each comprises 9 nodes; in 3D, the SUDS 
molecule increases to 27 nodes, whereas the QUDS size is only 13 nodes. Thus 
QUDS possesses a nearby two-fold advantage in coefficient assembly time and 
storage requirements. 

Also, it should be noted that Bounding Scheme 1 is inapplicable to the SUDS in 
3D, due to the fact that the coefficients of the non-principal, or "corner" 
nodes, may become negative. (Bounding Scheme 2 is of course still applicable). 

Relatively few applications have been reported using SUDS and QUDS in flow 
field calculations. Fewer still have been made using BSUDS1 and none using 
BSUDS2 or BQUDS2. 
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The studies of greatest relevance here are those of Castro (Ref. 25 ), 
Leschzeiner and Rodi (Ref. 26 ), Han, Launder and Humphrey (Ref. 27 ), 
Leschzeiner and Launder (Ref. 28 ) and Lai (Ref. 23 ). Their findings may be 

summarized as follows: 

i) No major problems were reported in solving the SUDS-based equations, but 
difficulties were encountered in obtaining solutions when QUDS was used. 
The exception is the recent work of (Ref. 27 ) where a particular 
arrangement of the coefficient matrix was found which produced stable ADI 
solutions in comparable numbers of iterations to those for the UDS 
calculations of the same problem. 

ii) The accuracies of the SUDS and QUDS results were similar in most 

instances, with the notable exception of the test problem of laminar flow 
in a square cavity with a sliding wall. Here SUDS was scarcely better 
than UDS, whereas QUDS was markedly superior to both. 

iii) BSUDS1 has been shown by Lai (Ref. 23 ) and Benadecker et al (Ref. 29 ) 
to perform satisfactorily as regards convergence and to be markedly more 
accurate than UDS (again with the exception noted in (i) above). 

The poor performance of SUDS referred to in paragraph ii) was attributed by 
Lai (Ref. 23 ) to the failure of this scheme to account properly for the 
effects of pressure gradients, which act like "sources" in the momentum 
equations and are particularly important in the cavity problem. This 
deficiency, which did not appear to manifest itself in other problems 
involving strong pressure gradients, arises because SUDS tries to preserve 
velocity, rather than total pressure, along streamlines. 

4.9 CONCLUSIONS 

1) Of the schemes examined, BSUDS2 and BQUDS2 appear to offer the most 
promise for 3D applications. 

2) The cost-effectiveness of QUDS and BQUDS2 is currently impaired by the 
lack of an efficient and reliable solver. Once a suitable solver is found, 
then the accuracy and compactness of the 3D computational molecule of 
these schemes is likely to render them superior to SUDS and its 
derivatives. 

3) The ADS, although of higher formal order than SUDS and QUDS, is 
non-conservative, and has been rejected for this reason. 

4) Spline schemes of the conventional centered variety like the CSS appear to 
be susceptible to severe overshoot problems and are also non-conservative. 
For these reasons they are not recommended at their present stage of 
development. Similar comments apply to the GRHS Hermitian scheme. 

5) The flux-blending strategies for eliminating overshoots seem ideally 
suited for implicit calculations and are likely to see further development 
and application. 
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4,10 SELECTION OF TWO DIFFERENCING SCHEMES 

No scheme seems to have emerged from this evaluation with a clear-cut margin 
of superiority over the others, SUDS offers a more favorable matrix structure 
from the solver point of view, especially when bounded as in BSUDS1/2. On the 
other hand, it has not performed well in isolated instances and it does 
generate a large computational molecule in 3D, 

QUDS has performed well in flow field applications as regards accuracy 
(although it too produced some isolated poor results in the present scalar 
transport tests). Further, it has the compelling attraction of a comparatively 
small 3D molecule. Although the basic scheme can produce significant 
overshoots, it seems amenable to bounding by the flux blending strategy 2. The 
main argument against QUDS is the lack of a wholly suitable solver. 

Both QUDS and SUDS, bounded and unbounded, show a flow angle dependence 
regarding their accuracy. QUDS seems to be more accurate at small angles 
whereas SUDS has a greater advantage at large angles. It should be noted that 
either of these two schemes will be an improvement over the hybrid 
differencing scheme at any angle. 

Thus, SUDS and QUDS were selected for evaluation in the two-dimensional 
version of TEACH. Although it was desirable to implement them with flux 
blending, it was not essential for testing these schemes in realistic flows as 
required in this contract. The unbounded version of QUDS was selected since it 
could be coded relatively quickly in TEACH and permit extensive testing of the 
method. BSUDS1 and BSUDS2 were also selected since they offered accuracy 
improvements comparable to those obtainable using QUDS, nonphysical 
oscillation could be avoided using flux-blending, and the computational 
molecule was reasonably compatible with the current solver. Based on the 
results of the test cases, it was anticipated that BSUDS2 would be more 
accurate than BSUDS1. It has been noted, however, that BSUDS1 is not suitable 
for use in the 3D-TEACH code. 
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5.0 DERIVATION OF SELECTED SCHEMES FOR TEACH 


In this section, the details of the derivations for the two improved 
finite-difference methods incorporated into the two-dimensional version of 
TEACH are presented. In Section 5.1, the derivation for the skew-upwind 
differencing scheme using the second bounding method (BSUDS2) is given. In 
Section 5.2, the derivation for the Quadratic Upwind Differencing Scheme 
(QUDS) is given. BSUDS2 was also selected for use in the three dimensional 
version of TEACH. The extension of the material in Section 5.1 for use in the 
3D-TEACH computer program is presented in Appendix B. 

5.1 THE BOUNDED SKEW-UPWIND DIFFERENCING SCHEME 

In this section, the bounded skew-upwind differencing scheme is described. 
First, a brief review of the flux form of the equations of motion is 
presented. Second, a detailed description of the finite-difference form of the 
flux contribution to a representative face of a typical control volume is 
given; the derivation of the flux contributions to the other faces is then 
outlined. Third, the resulting coefficients for the finite-difference 
equations representing the total flux (and sources) are presented. Fourth, the 
results of applying the boundary conditions are shown in a manner consistent 
with the foregoing flux representation. Fifth, the bounding scheme for the 
coefficients is detailed. 

5.1.1 Flux Form of the Equations of Motion 

The equations of motion for both laminar flow and (time-averaged) turbulent 
flow can be written in similar fashion for all of the dependent variables: 


~ — (pu$) + *— t - — (r^pv<j>) 
9x r° 3r 



7 * 


hi 


v£)* 


(5.1) 


where 6=0 for two-dimensional (planar) flow and 6=1 for axisymmetric flow. 
The variable i represents any of the dependent variables (e.g., the velocity 
components u, v, w, mixture fraction, turbulence kinetic energy and turbulence 
energy dissipation). The exchange coefficient, rtfs represents the sum of 
both laminar and turbulent contributions and is interpreted as the effective 
viscosity for i = u, v, w, the effective diffusivity for d = mixture fraction, 
etc. Sj} is a generalized source term. 

Equation (5.1) is integrated over a control volume appropriate for each 
dependent variable i and, after some manipulation, the finite-difference 
equivalent form of Equation (5.1) is obtained. The control volumes are defined 
using an orthogonal grid formed by the intersection of coordinate lines in 
both the axial and radial coordinate direction. For planar (two-dimensional) 
problems, a unit thickness and uniform properties in the direction normal to 
the plane of the grid are assumed. The intersection of the grid lines form the 
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grid nodes at which all flow properties except the axial (u) and radial (v) 
velocities are calculated; i.e., all scalars and the tangential velocity 
component (w) for swirling flows. The axial velocity is calculated using a 
second grid with grid nodes located midway between the scalar grid nodes in 
the axial direction and coincident with the scalar grid nodes in radial 
position. The radial velocity is calculated using a third grid with grid nodes 
located midway between the scalar grid nodes in the radial direction and 
coincident with the scalar grid nodes in axial position. Directions in the 
grids are identified as north, south, east and west. The grid system for the 
scalars is shown in Figure 5-1. The positions of the u and v storage locations 
(grid systems) relative to the scalar grid are shown in Figure 5-2. 



Figure 5-1 Grid System for Scalars 


1-2 1-1 I •+ 1 l + 2 
















-H 

^ U VI 

U U. J) 

f 

ELOCITY CONTROL 

VOLUME 

k — 




t $ < 

1 

1 

— M 
1 
1 

L_ A 

v- — i ' 

1 

1 vv 

,vu.ji \/ 

i 

i 

_j 

'ELOCITY CONTROI 

L VOLUME 



>» — — f 

► 1 ■■ 










Figure 5-2 Storage Locations for Axial (U) and Radial (V) Velocities 
Relative to Scalar ( i ) Grid System 
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Control volumes for each scalar and the tangential velocity are defined by 
planes located midway between the scalar grid nodes, as shown in Figure 5-3. 
Normal velocity components, therefore, lie along the boundaries of the control 
volumes. Generally, boundaries for the u or v control volumes include scalar 
grid-lines and are not necessarily located midway between velocity grid lines. 
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Figure 5-3 Scalar Control Volume 


The finite-difference form of Equation (5.1) is derived by integrating this 
equation over the appropriate control volume. For the scalars and the 
tangential velocity component, the grid systems defined above provide some 
computational convenience. Since the u and v velocities are stored midway 
between the scalar grid nodes, the convective fluxes from the north, south, 
east and west faces can be calculated without recourse to averaging any of 
these velocities. Also, the pressure gradient driving the flow can be computed 
without averaging pressures. 

In performing the integration over the control volume for each term in 
Equation (5.1), the mean- value theorem is employed and the source term is 
linearized in the vicinity of the center of the control volume (point P). 

After some manipulation, the finite-difference form of Equation (5.1) is 
obtained. 
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HU ‘ V» * Vn - C S*» " <*E ' V ‘ ®W ( *P ‘ V (S 2) 

- d b ( *n ‘ V ” D s ( *p " *s’ * <s » * s pV 

where Cg, C^, etc. are "convective coefficients" as defined below 
c £ - (pu) e a £ 

°N * (pv) n a n 

(5.3) 

C W “ (pu) v 8 w 
C s « (pv ) 6 a £ 

and a n = a s , a w = a e are the areas of the faces of the control volume. 


The "diffusion" coefficients are given by 



and ( ax ) e is the distance between points P and E, etc. (e.g., see Figure 
5-3). 

It is important to note that Equation (5.2) applies to all of the dependent 
variables although the appropriate grid must be used in each case to define 
the geometric parameters used in the calculation. Also, Equation (5.2) applies 
to all of the difference procedures considered in the program since each 
scheme is simply an alternative method for defining (interpolating for) the 
dependent variable at the faces of the control volume (e.g,, 0 e , jL, 0 n , 

$ s ). However, the diffusion terms are always represented by central 
differences. 
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It is convenient to define a total flux for each face of the control volume as 
the sum of a convective flux and a diffusive flux such that: 


where 


F e “ ^ e " d e ~ 

F w * Mw " ®w (lJ, p " ^ 

F n * ^ 

F s * C s*e “ D S ( *P " V 


(5.5) 


(5.6) 


In the following section, the skew differencing procedure will be used to 
calculate the values of the dependent variables at the faces of the control 
volume. As a result. Equation (5.5) will include not only the values of t at 
the "normal." or main, qrid node locations (E„ W. N, S and P) but also at the 
corner locations (NE, SE, MW, and SW). The finite-difference form for Equation 
(5.5) will be shown to be: 


VP • V'E + Vw + VN + V S + ^NE + A SE*SE 
4 Vnw 4 A S w*sv 4 s u 4 s p*p 


(5.7) 


5.1.2 Calculation of the Fluxes 

Recall that the equation of motion. Equation (5.1), can be written in terms of 
fluxes to each face of the control volume. Equation (5.5). In this section, a 
procedure will be described to calculate fluxes, F e , F w , F n , and F s . 

The derivation of F w , the flux to the west face of a typical scalar control 
volume, is given in detail. The derivation of the other fluxes and of the 
fluxes for the u and v velocity components are outlined. 

Consider the control volume shown in Figure 5-4. At present, it is assumed 
that the velocity vector is oriented as shown; i.e., the u and v components 
are non-negative. The flux to the west face is given by: 

F v = - D w ($ p - <f> w ) (5.8) 
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Figure 5-4 Control Volume for Skew-Upwind Differencing for West Face of 
Control Volume and Positive Velocity Components 


For central -differencing (CD), the value of the dependent variable at the west 
face, tf w , is given by linear interpolation between 6 at the W and P grid 
nodes. 

♦w = (1 " V* *W + °V*P (5<9) 

where the interpolating factor is 


a 


w 


X w - Xy 
Xp " 


(5.10) 


In the 2D-TEACH computer program, a = 0.5 for each face of the control volume 
for each scalar variable and the tangential velocity component since the 
control volume faces are located midway between scalar grid nodes. For the 
axial and radial velocity components, however, the a for each face may assume 
other values. The central difference form of the flux at the west face is then 


f«CD 



4 


♦w * ( ~ 1 > * P 


(5.11) 


where the Peclet number at this face is given by: 



C W /D W 


(5.12) 
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The upwind difference [UD) form for the flux, at the west face is obtained by 
setting a w to zero and neglecting the diffusion term: 


Fy i m 


*v 


(5.12) 


Equation (S T 8) Is also the starting point for the stew upwind differencing 
(SUD) scheme. The value of 6 at the west face of the control volume t is 
determined by extrapolating the velocity vector upstream to the point w r which 
lies along the grid line connecting the west and southwest nodes (see Figure 
5-4) to give" 


e w ' - U - K*} + 4 sw (5,14) 

where the skew Interpolation factor is the ratio of the vertical rlf stance 
between w 1 and SV to the vertical distance between W and SU: 


. 1 ok 

K ’* ^ 


IS. 15) 


For very large flow angles (skewing) relative to the coordinate directions, 
t w will exceed unity end w J will be defined in terms of 6 at the SW and 5 
nodes, however it is known that this approach can yield negative coefficients 
at the corner nodes [NW, SU, HE* SE} which can in turn produce oscillations fn 
the solution. To assure that the coefficients for the corner nodes are 
non -negative, them 


^ - max ^ 


The use of absolute values in Equation (5*16} permits this equation to be used 
to define k w for all velocity components at the west face. 

In his original development of the skew upwind differencing approximation, 
Rafthby [Ref, 11} assumed that # w ■= Thus* 


tx \ 
a., i 


rs.is} 





(5.17) 
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It is desirable to use tine central difference procedure for small values of 
the grid Peclet number and the skew upwind differencing method for large 
values of the grid Peclet number. It is also desirable that those two 
formulations produce a continuous transition at the transition Peclet number* 
which in the present case 1st 



for the scalar grid System used In the 20-TEACII computer program, Po* * 2 h At 
the transition Peclet number, the central difference result (Equation 5.11) is: 


F *ce _ i |w 

while the skew upwind differencing method (Equation 5 f 17) fields: 


1^2 * hL , t * p - 


(5.50) 


Noting the definition given by Equation (5.14], It is clear that these two 
results are not equal. 

The fluxes at the transition Peclet number can be rcade equal by noting 
(contrary to the assumption made by Raithby) that d w and i w ' are related 
byi 





uE * 


so that Equation (5.17) becomes 


(5.21) 


^51JP 
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Writing the central difference result In terms of the flux definition. 
Equation 5.8, then: 



P «V*V 



(5.23) 


Clearly* these two fluxes will be equal at the transition Peclet number If a 
correction, Pe£ [ S^AS* Is added to the skew upwfnd differencing 
flux, Equation 5.17, to obtain; 


°u 



The derivative [asi/9s) w can be computed by; 


15 . 24 ) 


£4 _ \ ~ **' [5.25) 

31 " A* 

Then, using Equation [5.&K [5,14) and (5,15), Equation (5.24) becorses 





k v ( *w “ *BW^ 


[ 5 , 26 ] 


At the transition Peclet number* the fluxes calculated by central differencing 
(Equation 5,19) and skew upwind differencing (Equation 5,26) are equal. 

It will be recalled that the above result for the skew upwind differencing 
flux at the west face of the control volume was derived for non -negative 
values of the axial and radial velocities. Similar results can be derived for 
other combinations of u and v by consistent application of the process leading 
to Equation (5,26). The result is & general expression for the flux as 
calculated using skew upwind differencing; 

“jp 2 - P e [o'* +U + fl - e“> * J 
r H \ v w w p 

- ««„ - v *.(■£{*« -°v ♦» - <» - <» v) 

‘"'OU-Vs-' 1 - O }) (J.ii 





The pan meters,. and are switch; that indicate the direct ion 
of the components of the local How velocities. 



fs.eai 


« v . i /i - t 

v 1 \ |Vj } 1 

Each of these parameters has a value of unity if the velocity component is 
positive (or, by convention, non-negative 1 and zero If it is negative. The 
transition Peclet number is now given by the general result. 


p* . i ts.m 

e V Ho ” 

°w * “ “v^ 

Ln the hybrid differencing procedure, the more accurate centra? differencing 
formulation (Equation 3,11) is used when the Peclet number Is less than the 
transition value while the less accurate ^ but stable* upwind differencing 
result (the generalization of Equation (5,13)) 


Du 


K % * U 


<> S 1 


(5.31) 


is used when the Peclet number is greater than the transition value. 

Originally, It was believed that a similar hybrid procedure could be developed 
for skew upwind differencing with Equation [ 5, II > used for Pe < Pe* and 
Equation (5,57) used for Pe > Pe*. However, this approach proved to be 
unworkable since some of the coefficients derived from this hybrid formulation 
for use in Equation (5.7) could be negative. As an alternative, a flux 
blending scheme (to be described In Section 5.1.5) is used In which a weighted 
average of the upwind differencing and the skew upwind differencing fluxes is 
used. The weighting (blending) factor, Tp . is chosen in such a way as to assure 
boundedness. The bounded skew- upwind differencing iBSUD) flux is defined by: 

(5.32] 


f *bSUD 

°W 




TwuB 
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with the lighting factor restricted to the range* 0 l a As a 

consequence of this def 1 nl tl on ' ~ ~ 


r|f B5LT> 

"to 


* % l % % 4 ° " 0 V 1 " <T *v " P 4 J ^ 

\°l {*W ~ ff w f SW ' ^ 

+ * d v +a - U **}] 


f 5.335 


finally, a bounded skew hybrid differencing (E$H&) formulation can be defined 
as: 


Fw asnE- 



[1 - V T 


U BSU1> 


(5,34) 


where \ is- permitted to assume only two values: 1 for central 

differencing (Pe < Pe*) and Q for bounded skew-upwind differencing 
(Pe > Pe*], Equations (5,33) and (5,34) are the basic working relationships 
used to determine the flu* at the west face of the control volume. The 
contributions of the flu* to each of the coefficients In Equation f 5. 7 ) can be 
tmnediately identified by using Equations {5.33] and (5.34) In Equation (5,5). 

Equations analogous to Equation (5.33) can be derived for the other three 
faces of the control volume in exactly the same manner as employed, above. 
However, the results can be obtained by inspection as follows. 


?l 



The east face flux fs obtained by translating the nodal subscripts eastward 
such that: 


West Subscript 


East Subscript 


becomes 

becomes 

becomes 

becomes 

becomes 

becomes 


Of course, the lower case subscript "V becomes "e, " 

The south face flux is obtained from tbe west flux (Equation S. 33 J by rotating 
the nodal subscripts counterclockwise through W degrees. 


West Subscript 


East Subscript 


V 

SK 

NW 

$ 

N 

u 

l - v u 

(TV 

1 - ff'* 


becomes 

hecomes 

becomes. 

becomes 

becomes 

becomes 

becomes 

becomes 

becomes 


£ 

^E 

SW 

F 

W 

<jV 

1 ^ ff V 
l . CT u 


The north face flux is derived by translating the south flux result northward, 


South Subscript 


Worth Subscript 


becomes 
bat rarces 
becomes 
becomes 
becomes 
becomes 


The results of the manipulations are summarized in Tables 5-1 through 5-111. 
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table: f-i 


DEFINITION OF OOtINDfD SKFW HYBRID DIFFERENCING FLUXES 
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TABLE 5- I I 


CENTRAL DIFFERENCING FLUKES 


vest 


CD 

»w 


(L -a u > * 1 ] , M . (v«v‘ ») ♦, 


iitE 


*CP 

u £ 


P„. (1 - O * 1 
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TABLE 5-1 II 


BOUNDED SKEW-UPWIND DIFFERENCING FLUKES 


% l % +V + a - <0 * P 1 ‘ < p e w - TJl, 

j 

“i *H- *SH~ (1 - <C> *™ 

'll 


P e c « E 1 - (T < 


*P ‘ °l *s - C > - 'I 5 *N 


* (t ‘ +E ‘ "I *SE ' (1 ' V 


-ff 5 * p =« K * * » - °:> v - <*«„ - v v. 
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% k - (l - O *. 


*' SE s T SH 


+ (I - ^ ~ (1 - i£) *e “ ^ J 


ff n *F " (1 " O *E - u“ #w 


+ n - I L - (1 - t? U > £ - cr v 1 1 1 

n N tv 7 W KE n Vrj 



5.1.3 Calculation the Coefficients for the Finite-difference Form of the 

Equations of Motion 

The finite-difference form of the equations of motion, fe.g., Equation (5+?)) 
can be derived directly from the flux Information presented in Tables 5-1 
through 5- III and the sign conventions determined frwn Equation f 5 .5} . The 
resulting expressions will contain the unknown blending factor, y r The 
blending strategy requires that the terms In the equations for the 
coefficients most responsible for producing negative coefficients he isolated 
so that appropriate values for T can be determined. Furthermore, the 
coefficients for the control volumes adjacent to the physical boundaries of 
the flow may require modification to incorporate the effect of the boundary 
conditions. Thus, to simplify manipulation and modiftcatl on* some additional 
notation will bo defined. 

Let the center of the control volume (point P} be located at the Ith axial 
position and Jth radial position. The flux contributions (the components of 
the total flux) to the east face are denoted as E?(1,J), E3(I,J) and 

the flux contributions to the north face are denoted as N1(1 T J) N£(U), 
N3(l t J). These flux contributions are defined as follows. 

Central Differencing 

ZUU) - Of - *pCe 

E2(I,J) = EIU,JI * C E 

E3( I , J) = 0 

HUM! - Djj - *r£u 

NZ(I,d) p HI [ I P J J 4 C^ 

- 0 

Bounded Skew-Upwind Differencing 


SHl,J> - to" - i> 
E2C5.J) - 6 K I , J) + c £ 

- k.tC, - p.’oj.) 

' 1) 

- HUlpJ) * C H 
F3( I , J) ^ kjtCn - 


(5.35) 

15.3C) 

(5*3?) 

(5.30) 
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The use of central versus bounded skew-upwind differencing is determined by 
the value of the Peclet number at each face. The parapieters Cp, Dr, a e , 

, .are local values; the subscripts fl*Jl nave been omitted in the 
interest of readability The corresponding flu* contributions at the west face 
are given immediately by E1{]-1,J), E3(I-1,J) and the flux 

contributions at the south face are 101(1**1-1), M3fl*J-l). 

The coefficients of the finite -difference form of the equations of [notion may 
then be defined in terms of these flu* contributions. The results are 
presented in Table 5-IV. The coefficient Ap can be shown to be equal to The 
sum of the other eight coefficients when use is made of the mass continuity 
restriction: 


*■ m 0 (5,39) 

The boundary conditions {Section 5.1.4) and the blending scheme {Section 
5.1.5) can be applied directly to the flu* contributions so that the results 
shown in Table 5- I v are general, 

5.1.4 Boundary Conditions 

The boundary conditions are applied to the flu* contributions as defined by 
Equations (5.35) through (5.38} SO that th& set Of coefficients (e.g., Table 
5-1 V) for each dependent variable f$ the same throughout the computational 
domain. This consistency simplifies the application of both (!) the algorithm 
for solving the set of simultaneous equations for each variable and (?) the 
bounding procedure (see Section 5.1.5), 

5.L.4J Scalars 

A typical control: volume for a scalar variable is shown in Figure 5-5 (which 
is based In part on Figure 5-1), It will be recalled that the axial velocities 
are stored midway between the scalar gridlines in the axial direction; 
therefore, the axial velocities are stored at the midpoints of the oast and 
west faces of the scalar control Volume, Similarly, the radial velocities are 
stored at the mid points of the north and sooth faces of the scalar control 
vol umes . 

The extreme physical boundaries of the flow are located midway between the 
most extreme gridline and its nearest neighboring gridline In each direction. 
The gridlines are numbered from 1^1 through I=NI in the axial direction and 
from 0-1 through d=NJ In the radial direction. Thus, the interior nodes are 
numbered from 1-2 through lll-l In the axial direction and from J=? through 
NJ-i In the radial direction, In the following discussion, the boundary 
conditions are applied to the various faces of the control volume for an 
Interior node {1*J), as shown in Figure 5-5, to produce modified east-west 
flu* contributions {El, E2* E3} or north- south flu* contributions (Ml, NJ\ N3) 
for use in calculating the coefficients of the finite-difference equations for 
the scalars (e.g,.. Table 5-1 V). In the process of introducing the effects of 
the boundary conditions for the Interior nodes. It will be seen that some of 
the flux contributions at adjacent nodes are also affected. 
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TABLE 


CQ- EFFICIENTS OF THE FINITE - D I fF£ REWCE EQUATIONS 


Aytt.J) = E2U-1.J) - Tr v E3 Jj_j _ j 

+ y b Hi “■ o v > c u Ni] T T-1 

p hd IjJ I 

M n l { 0^ H3] x j 
■ii n n 1 1 J 


A E a,J) = E 1 < 1 „J> + y e l(l - v") E 3 1 T + T r<] - cr v Hi - W 3 ], t . 

■“ ^ G IjJfi ft fct 

-Y_lc* (1 " O J 

cd n I , J 

AgU.J? = N 2 (r h J-l) - V [a v m, . , * T w [Cl - o U ) o V 133 . , 

° “ n * « I-l,J 

■ % K % “ J i., 


a^i.JI - nKi , j!) * 1(1 - g v > Nil + [Cl - P u )(t - e v > ESI t _ 

11 w n IjJ c e I-ljJ 

- 1. 1 .; a - •:» **!„ 


A s w d t J) 


"'mj * t - 


[o“ o V »] ... 

n n 1 P J“* 


AgjClpJ) - > t? V U - <■“? M 3 L , , - V e [<l - " U > cr V E 3 ] 

s n n - e e £,J 

AyatlpJ) - y [o’ J (3 ^ o V ) E 3 ] " T n K 1 " ^ 3 ** 1 , t 

1 “ e e 1-1 ,J r. n 

AmClpJ) * -T 4 tfl “ - S V ) E33 - Y n ttl - c V H\ - o y >]N3 

“ e e. t T,J ” n n I f J 


NOTE: The blending factors t 6J v s and arE the factors 

appropriate for the (IpJ) node as determined by the procedure discussed 
In Section 5,1'$. 
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There are Four types of boundary to be considered: 
fl) specified wall, 

(Z) specified inlet, 

(3} ex is of symmetry* 

(4} unspecified opening, 

The application of these boundary conditions to modify the flux contributions 
used to determine the finite-difference coefficients for a node (1-1*3) having 
a physical boundary on the west face of Its control volume will be described 
in detail, The modifications to the coefficients due to the other physical 
boundaries are similar and will be Summer f ted. 

n est Boundary 

for the specified wall boundary condition* the velocity normal to the wall at 
the west boundary (the west face of the control volume for the node [I.J)) Is 
zeroi velocity components In the radial or tangential direction may be nonzero 
(e.g. , a moving wall]. The zero normal (ti) velocity, which for this case Is 
Stored at the axial location of the west boundary* results fn the use of the 
central difference formulae for the flux contributions, El, F? and £3, for the 
west face of the node U 3 lT) + These flux contributions are computed and stored 
as east face contributions for the node see Equation f5,35]. 

Specifically* for g specified wall located at the west face of the control 
volume for the node 


E1(I-1*J) = DE( I-L,J) 

F2(I-1*J) * Eld-l.J) (5.40) 

E3fl-l,d) * 0 

for a specified wall, the value of the scalar at the wall may be either known 
(e.g.* turbulence kinetic energy = 0) or unknown (e + g,, temperature). Since 
the variation of the scalar from its wall value to the value at the node (I ,J) 
is generally non-linear* it is prudent In every case to decouple the wall 
values from the system of equations. The influence of the wall can be 
reflected by using appropriate wall functions to compute the source terms for 
the scalars. From Table 5-IY* it can be seen that the known or unknown wall 
values of the scalar can be excluded If the coefficients A w . A 5W . A^y 
are zero for the node (1,0). These coefficients depend upon the values of: 

a) n (I-M) 

b) £3 (I-I*J) 

c) M3 U,J) - H3 at the north face of the control voluete for node 
(I*J) 

d) N3 - N3 at the south face of the control volume for node 

(M) 

Inspection of the coefficients listed in Table 5-IY indicates that N3 ( I , J ) Is 
excluded automatically unless the axial and radial velocities are positive at 
the north face* and N3(I,J-1) is excluded unless the axial velocity is 
positive and the radial velocity is negative at the south face of the control 
volume for the node (I*J). 
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The necessary modifications can be summarized as follows. 


El(I-LJ) = o 
N3(1,J) = D 
H3(l,J-n = Q 


(5.41) 


It follows Immediately from Equation (5. 40) that: 

£2(1-1, J) = 0 

E3(I-1 ,J) = Q (5.4?) 

and therefore p Ay, Ajy, and Ayy will be zero for the node U,JL 

For a specified opening, the value of the scalar is known at the west face of 
the control volume for the node In general, the cell Peclet number is 

greater than the transition value whenever the boundary represents a specified 
opening. Therefore* the flu* contributions EL 11 and E3 at the node (I-LJ) 
are computed using the skew upwind differencing formulae. Equation ( 5*37) * 
however* since the boundary is coincident with the west face of the control 
volume for the node there Is no skewing of the flow at this face and, 

therefore, E3(I-1,JJ is set to zero. 

The determination of the coefficient Ay for the node (I,J) also requires the 
computation of the north face and south face ft us contributions, N3[I,J) and 
respectively. Since the specified opening is coincident with the 
west face of the control volume for the node (I*J) and since this face is 
located midway between the nodes ( I -1* J ) and it Is necessary to double 

the value of the skewing Interpolation factors used to compute M3(I,J) and 
N3([.J-1); l.e.. 


k n = max 

k s b max (?*k s , 


1 . 0 ) 

LOi 


(5.43) 


where the values of k^ and k s on the right-hand side of Equation (5.43) 
are calculated in accordance with the normal distance between gridlines 1-1 
and L 


If the west fate ef the control volume for the node [l,J] Is. an axis of 
symmetry > then the flu* normal to this face Is zero. By definition, the 
unknown values of the scalar at the nodes (I-Ld) and (I,J) are equal;, the 
gradient of the scalar normal to the axis vanishes. Thus, an axis of symmetry 
Is mathematically Identical to a specified wall boundary condition with zero 
gradient (source) at the wall.. The modifications bo the flux contributions are 
given by Equations [5,42) and [5,43). 

For an unspecified opening, ft is assumed that ( T) the flow is exiting through 
the opening* (2) the streamlines are parallel fh the vicinity cf the opening, 
and (3) the cell Peclet number exceeds the transition value. Determination of 
the coefficient Ay, Ajy, and Aiuy for the node (Ld) requires the 
calculation of the flux contributions E2U-LJL E3(I-1,JL H3(LJ) and 
K3( I „J-I), By the third assumption, the skew upwind differencing formul ae. 


Si 



Equations (5*37) and (5.38), are used, By the second assumption, E3(I-I*d), 

W3C I P J) and II3(J,J-1) are zero. By the first assumption, the flow direction 
switch 0 “ at the node U-1,J) Is zero so that is zero. Thus these 

assumptions are equivalent to: 

£S(M.J) = o 
E3(I-1*J) - 0 

U3(I,J) = 0 (5.44| 

N3(I»iM) * Q 

and it fellows that Ay, Agy* and Ayy for the node (1,J] vanish. Thus, 
the unknown values of the scalars at the unspecified opening are decoupled 
from the system of equations. 

It Is seen that the modifications to the flux contributions for the specified 
wall, axis of symmetry, and unspecified opening boundary conditions are 
Identical and* therefore, the coefficients Ay, Ajy, and A(jy for the node 
{1,J) vanish for these boundary conditions. 

East Boundary 

If the east face of the control volume for the node (I*Jl Is a physical 
boundary, it Is necessary to modify the flux contributions Elfl.J)* E31I.3K 
W3C 1,0) and tJ3( used to determine the coefficients Af, A^r, and 

A|je for the node H.J). for reasons similar to those given In the discussion 
of the west boundary. In the case of the specified wall* axis of symmetry, or 
unspecified opening boundary conditions* It is necessary that: 

Ein*jj = o 

E3(I,J) = 0 (5.45) 

N3U*Jl = 0 
N3(I,J-1) - 0 

so that %, Asf, and % for the node (1,3] are zero. For the specified 
opening, it Is sufficient that: 

E3U,J) ^ 0 

k n • max (Z*fc n * 1.0) (5,45) 

lt s . max (2*k SJ 1.0) 

The last two conditions guarantee that the flux contributions N3(I*J) and 
H3f I P J-1 ) are calculated correctly. 



South' B&unflary 

If the south fat& Of the control volume for the node (1,J) Is a specified 
wall, an axis of syumietry, or an mnspeclf led opening, then it is necessary 
that: 


N3(l ,<M) = 0 

E3U-1.J) =* D (5. *7] 

£ 3 ( 1 , 3 ) - 0 

so that the coefficients Aj., Agy a and Agg 'finish. For a specified 
opening, It is necessary that: 

W3U.J-1) t 0 

k w = mas (2*k*, 1.0} (5.46) 

kp = mas (2*k e , 1.0] 

The last two conditions assure that and £3 f I * J > are computed 

properly. 

Worth boundary 

If the north face of the control volume for the node (1,0) represents a 
specified wall, an axis of symmetry, or an unspecified opening, then ft is 
necessary that; 

HI (1*0) = 0 
N3(I,J) = 0 

E3(M.J) p 0 (5,49) 

E3fl,j| = 0 

so that the coefficients Ay. A^y. and A^e vanish. For a specified 
opening, it is necessary that: 

>13(1,0] o o 

k w = max [Z*fc w , 1.0} (5.50) 

k e - max (2*Jt e , 1.0} 

The last two conditions guarantee that E3(I-1,J) and £3(1,3) are calculated 
correctly. 

5. 1,4.2 Axial Velocity 

A typical control volume for the axial velocity, u, Is shown in Figure 5-6, 

The indices I and J refer to the scalar grid system. The index I 1 refers to 
the axial velocity grid system; note that this Index varies between I'«2 and 
Nl. The east and west faces of the axial velocity control volumes are 
coincident with the vertical scalar gridlines. The axial storage locations of 
physical boundaries are the same as the axial locations of the axial 
velocities; the Inverse is, of course, not true. The radial locations of the u 
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velocity and scalar control volumes are identical- The radial velocity, v. fs 
stored (effectively) at the corners of the w -velocity control volume. The 
system of simultaneous equations for the axial velocity field describes the 
4X1 &1 velocity distribution from 1 1 to Nl-1 and J=£ to MJ-l. However, the 
coefficients for thfs system Of equations require the determination of some 
flux contributions along the gridlines l'=2 to I 'nNl-1 and 0*1 to 0=MJ-1. 


(D sk^gf like roh jttskal proum 



Figure 5-6 Scalar and d-Velocity Grid Systems Showing U- Velocity Control 
Volume Near Physical Boundaries 
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Ne&t Boundary 


If * physical bounds ry f$ located to the west of the axial Velocity control 
volume For the node (T*J)* then the physical hound* ry is located along the 
gridline, I'-l. if the physio*! boundary represents either a specified wall or 
an axis of symmetry, then is zerOJ, if it represents a specified 

opening* then u(I , -l*J} Is a known, specified value. If the physical boundary 
represents an unspecified opening* then it fs assumed that the velocity 
utl'-ijJ) fs known from the previous Iteration using a special procedure 
described below. In every case, it is seen that the axial velocity is known at 
the node { I J - 1 , J ) + Therefore, the axial velocity at the west face of the 
control volume for the node can be computed in the same manner as fs 

used for any axial velocity interior node. 

The radial velocities for the west face of the control volume are stored at 
the corners of the axial Velocity control volume (see Figure £-6), Therefore 
the radial velocity at the center of the west face can be interpolated 
directly from these data. Thus, one can obtain the values of u and v and the 
flux contributions can be computed using Equation (£.35) or (5*371 in the 
usual manner. Thus, it is seen that the application of the boundary conditions 
at the node {L,J} requires no modification to the computation of EJ , £2 and E3 
provided that the correct axial velocities are stored at the node (I'-I^J). 

In the discussion on scalar variables, it was stated that, at an unspecified 
opening, it is assumed that (II the flow is exiting through the opening* (21 
the streamlines are parallel in the vicinity of the opening, and (3) the cell 
Peolet number exceeds the transition value. These assumptions are also made 
for the determination of the axial velocity at an unspecified opening. It is 
possible to estimate the axial velocity at the node using these 

assumptions together with the mass conservation restriction for the total mass 
flow rate and a relationship for the axial velocity distribution at an 
unspecified opening based upon the second assumption and conservation of 
radial momentum. Since this procedure is not peculiar to the bounded skew 
hybrid differencing procedure but is part of the general procedures used for 
all finite -difference schemes in the TEACH computer program, it will not he 
described in this section. The procedure produces an estimate for the axial 
velocity at the node (I , -1,J) for the unspecified opening boundary condition 
using the solution obtained from the previous iteration. 

East Boundary 


If a physical houndary is located Immediately to the east of the axial 
velocity control volume for the node (I’J), then the physical boundary Is 
located along the gridline, r*l. The various boundary conditions affect the 
value of the axial velocity u[I'*l,J) in a manner similar to that described 
for the west boundary. 
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South Boundary 

The redial positions of the scalar and axial velocity grid systems are the 
same, Therefore, If the south face of the axial velocity contra! volume for 
the node (I'jJ) is also a physical boundary , then the modifications of the 
flux contributions are similar to those made far the scalar with a physical 
boundary along its south face. Equation (5. *7) or (5.4B), In the present esse, 
the storage location for the radial position of the axial velocity is always 
coincident with the radial position of the axial velocity controT volume. 
Therefore^ the stewing interpolation factors, k w and k e , are always 
increased in accordance with Equation f 5*40) . 

If the physical (boundary at the south face of the axial velocity control 
volume represents a corner, then in accordance with the definition of the 
location of the physical boundaries relative to the scalar grid system, the 
corner extends to the midpoint of the south face. Therefore, the axial ftl s ) 
and radial iV s ) velocities used at this face to compute the flux 
contributions are determined from, a weighted average of velocities for the 
specified wall and an appropriate flow boundary condition. 

Horth Boundary 

If a physical boundary is located along the north face of the axial velocity 
control volume for the node ( 1 1 „ J ) , then the boundary conditions are used to 
modify the Flux contributions in a manner similar to that used for the south 
boundary. The treatment for a corner is also similar. 

5.1. 4.3 Radial Velocity 

A typical control volume for the radial velocity, v* in the vicinity of the 
physical boundaries is shown In Figure 5-7, The indices I and J refer to the 
scalar grid system. The Index J r refers to the radial velocity system; note 
that this index varies between J'.Z and MJ. The axial location of the v and 
scalar control volumes are Identical. The axial velocity is stored 
{effectively) at the comers of the v -velocity control volume. The system of 
equations for the radial velocity field describes the radial velocity 
distribution for 1=2 to Nl-J and J'*3 to NJ-1, However, the coefficients for 
the system of equations require the detennl natl on of some flux contributions 
along the gridlines 1-2 to W3-1 and J'=2 to MJ - 1 , 

Because of the synptetry of the notation describing the axial and radial 
velocity distribution, the modifications to the parameters used to compute the 
flux contribution coefficients of the finite-difference equations for the 
radial velocity follow immediately from the discussion of the scalar and axial 
velocity boundary conditions, 

5. 1.4. 4 Tangential Velocity 

The axfal and radial locations of the gridlines for the tangential velocity. 
w > and the scalars are identical. It then follows that the boundary conditions 
are treated In a similar manner so the modifications discussed earlier in the 
case of the scalars apply to the tangential velocity as well. 
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Figure b-7 Scalar and V-Velocity Grid Systems Showing V-Yelodty Control 
Volume Hear Physical Boundaries 


5.1.5 The Bounding Scheme 

The calculation of the hounded skew- upwind differencing fluxes and, therefore, 
the determination of the coefficients to the finite-difference form of the 
equations of motion,, require that the hi ending factor, T , be determined. The 
blending factor specifies the respective contributions of the flu* computed 
using skew and upwind differencing. For example* at the west face of the 
typical control volueie,, Equation {5.32) states; 




’Arid * {1 ' 


{S.51] 


The coefficients Including the local blending factor are listed In Table S- IV,. 

It Is possible to show that the corner coefficients (A^y* A^f* A^y* 

AffE) are unconditionally non- negative. For example, consider the coefficient 

Hv* 


a £W ei,j> * t w * T - t v? 33 i,j-i 


(5 + 52) 
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Both -ty and tg are restricted to the range 0 _< t <. !■ if both the axial 
and radial velocities at the west face of thE control volume are positive, 
then both and o]( at this location are unity: In any other case, 
either Ur both) are zero. For positive u. and v velocities, the flu* 
contribution £3(l-t,J) at the west face is positive. Therefore* the west face 
flu* contribution to Ajy is always non-negative. The same reasoning* when 
applied to the south face, leads to 4 similar conclusion* Therefore* A™ Is 
unconditionally non- negative* The other corner coefficients are treated In the 
same manner. 

How, consider the four main coefficients (Ay, Af, A s , A N ). Frnm Table 
5-1 V t 


AwU.Jl * £3(1-1, J) - v w [0^E3] Jwl J 


Y s 1(1 -v V J o“N3] 

6 fl TV 1.J-] 


- Y b [tj'n'Vj) 

n n I,J 


f 5.53) 


From the definitions of the flux contributions, Equations (5.35) through 
(5.3B)* it Is evident that; 

(1) Is positive if the axial velocity Is positive and It Is 

zero otherwise: 

(£) 0 ^ and E3 at the west face (I-l.J] are positive if the axial 

velocity Is positive but is zero if u is negative; 

(3) from the definitions of E2 and £3 and the fact that < -> it is 
therefore concluded that the first two terms In Equation f 5 , 53) must 
yield a non-negative result: 

(4) the third term is zero unless 0 Is positive and v is negative at the 
south face, in which case, it is negative; 

(5) the fourth term Is non- positive because the term within the brackets 
is always non-negative. 


Therefore, it is possible that the coefficient Ay Is negative. 



It is desirable that all of the coefficients of the finite-difference form of 
the equations of motion (Equation (5.7)) be non- negative for in this case the 
Value of the dependent variable d at the node P Is simply a weighted average 
of the values of d at the surrounding nodes exclusive of the somewhat 
complicating effects of local sources. A bounding scheme is a procedure to 
limit the values of the coefficients of the finite difference equations in 
such a manner as to produce this physically realistic result. Its principal 
computational advantage is to exclude under- and overshoots of the solution 
during the iterative procedure. These oscillations can produce severe 
numerical instability. 

The bounding procedure used herein is based upon the following sequence, 

(1) For each Iteration, the solution for the distribution Of the dependent 
variable fs Is obtained using the blending factors t determined during the 
previous Iteration; the blending factors are set to unity for the first 
Iteration, 

(?) For each point P In the variable field, the maximum and minimum values for 
$ are determined from the neighbors; 


W™ ( V V *£■ V V V V »SE) ,5 - m ’ 

( V V V V »NE' V *5E) I5 - 551 


It should be noted that the effects of sources are net included explicitly 
In calculating or since the explicit result is difficult to 
perform while avoiding double counting of the source effect. However, the 
effect of sources on the neighbors fs accounted for implicitly because the 
difference equations being solved Include these sources. In any case, the 
results obtained will be no worse than to force the use of upwind 
differencing by making ^ * 0 at the node where a peak in the profile of 
the variable being considered Is occurring. 

It should also he no tod. that one or more si may be excluded from the 
determination of or when the control volume under 
consideration Is adjacent to certain types of physical boundaries. For 
example* IF the x-axis is an axis of synstietry* then rfg Is excluded since 
It Is Identical to s) p and Its use could yield an inappropriate range of 
permissible values. 

(3) If ^Hiin < rfp < tfftiax* then the lota! value of the blending factor fs 
unaltered from Its current value, However, if dp Is outside of this 
range* then a new value of r must be determined. For this purpose each of 
the coefficients in Table 5-TV is written In the form: 


\ * \ 4 T V k 


ff, 5, HE, SE , » ^ 


iSM) 
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where It has bee* assumed that the blending factors for each of the faces 
of the local control volume are equal* Then Equation [&«7) -can be written 
as; 


II 4 ■■ 

* PVYt^i - y„ * E Vk 4 T 1 Vk +&M 


(6.S7) 


Now, If fr < rfmtr, then from this equation, 


Y > 


1 Mk ~ ^mCn E Aj, * Su * BpCain 
■^min * A k “ 1 A k 4 k 


ss.ss) 


tut if ip > imaxi then 


< LH** - W i»kt gfl 4 Sj t 


u yiMii 


W 1 *K ~ E Vk 


(5.591 


Of course, in practice the value of t is determined by use of the 
appropriate equality for Equation (5.5B) or Equation [E.5£). 

U) In this fashion, the local value of vff,J) is determined, The coefficients 
in Table 5- IV are then recomputed using the blending factors and the now 
distribution of *i is determined* For each face of the control volume at 
the point there are two values of t as determined for the two 

control volumes sharing the common face. Thus, a sufficient condition for 
assigning local values of the blending factors Is given by' 


\ * *i* U y (I T J)1 

'y “ It y 


* min [t ClpJ-lJ, Y 
Y - min fy Y C3,J+J>] 


(5,60) 


min 



This- procedure for determining the blending factors effectively limits the 
value of the dependent variable range to values of Its principal neighbors. 
However* It is still possible that some of the rrrflfn coefficients will be 
negative. More restrictive schemes tan be formulated fe,g + * blending strategy I) 
that guarantee that all coefficients of the finite -difference equation’s are 
non- negative but these procedures achieve this result by reducing the blending 
factor toward zero- the greater relative contribution of upwind differencing 
in this case produces larger amounts of numerical diffusion. Thus* the 
selection of a hi ending scheme represents a comer cm. iso between the undesirable 
effects of numerical instability and numerical diffusion. 


5.2 THE CHJ0S DIFFERENCING SCHEME 

In this section* the QUDS (for Quadratic Upwind Differencing Scheme) is 
described. This finite differencing method was developed by Leonard (Ref. 30 
and 31) and is based upon interpolating for the value of the dependent 
variable at each face of the control volume by using a second -degree 
(quadratic! polynomial biased toward the upwind direction as discussed below. 
The interpolated value is used to calculate the convective term in the 
governing equation for the dependent variable while central differencing is 
used to approximate the diffusion term. No bounding procedure was used with 
QUDS in the current program. 


5,2.1 Flux Form of the Equations of fiction 


It will be recalled from the discussion of the skew-upwind differencing scheme 
that the flux form of the governing equation is: 


F e - + F r _ - F g = Sy + Vp 

where the fluxes are calculated using; 
F e = c E “ &£((*£ " V 

F w “ r W p w _ % Ci p “ 

F n “ - V*N “ V 


f 5 , 61 S 


(5,62) 


F s = C S*s “ D S^ P * 

Then using quadratic interpolation for dp* tf w . and referring to 

Figure 5-8* it can be shown that the f imte-di fference form of Equation (5.61) 
is: 

Vfc = *E*E + + + *S e S 4 *EE*EE + VlArw + 15.63) 
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5. l.l Calculation of the Fluxes 

In this section, the calculation of the values of the dependent variable at 
the faces of the control volume shown in Figure 5-8 is described. The 
derivation of Fy, the flux at the west face of a typical scalar control 
volume, Is given in detail. The derivation of the other fluxes Is outlined. 

Recall that the flux to the west face is given byt 

Fy = (5.64) 

In its original form, QUDS required the use of a qeneral second-degree 
polynomial for the calculation of {see Ref. 31]; 

f w “ c O + c i* + c 2 x2 + C 3V + + c 5 J< y (5,65) 

However, It is argued in Ref. 31 that it is generally unnecessary to consider 
interpolation In the transverse direction so that for 0 w [and, similarly 
for d e ); 




Cg + Cj_X + CjX 2 


while for or d 5 : 

* = + c iV + c 2 y 2 


{ 5.653 


[5,67) 
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Figure 5-B Control Volume for the QUDS Scheme 
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'n calculating d w ,, the mathematically equivalent but more convenient 
quadratic form is used. 


t = o,, ■+ Ci lx - >tpJ + c 2 <k - ltpMs - 


rs.&ej 


Consider the situation for positive axial velocity u in which case the values 
of d at the grid nodes P* K and VH are used to evaluate cq* eg, and cj; 
that Is, the calculation of d w * designated d* In this case* Is biased 
toward the upwind node at X^. Then: 



°w ^»p . j *wtr* P (X'J-Xy) 

V^p (Xfcr ^ 1 £r*p " *w vr*w 


(5 r £9) 


For the case when the axial velocity Is negative, the coefficients in Equation 
(5.6S) are calculated using the values of 4 at the grid nodes W f P and E (that 
is, biased upwind from the west face) so that: 



V^n 

G + - — (Xy-Xj 

P V 5 ^ ^ 


■DEr-tp 


:_ e _ 

^ET^p 


V»p iyV <VV 

V^p ’ V*w 


15*70) 


Defining a flow velocity direction switch in the same manner as that used in 
the bounded skew- upwind differencing procedure, the value of 4 at the west 
face becomes: 


*w - 


S < + (1 - <£> +; 


(5.71) 


The value of d at the other faces of the control volume Is calculated 
analogously. For the east face, the values of d at the grid nodes W, P end E 
are used when the axial velocity Is positive and the values at the grid nodes 
P, £ and EE are used when the axial velocity Is negative. Since the 
interpolation for the south and north face Is formally Identical to that for 
the west and east face, respectively, the expression for jjl s (or dfi) can he 
Obtained directly from d« (or dg], 

5.2.3 Calculation of the Coefficients for the Finite-Difference Form of the 
Equations of Motion 

The finite -difference form of the equations of motion fe.g., Equation (5. 63)} 
can be obtained; by Inserting the expressions for jj w , ris, and d r 
Into Equation (5.62) and applying Equation (5. 61}, From equations like (£,69) 
and (5.70), It con be seen that all (J| T 1 £ ?, occur as the difference 
between 4-\ and dp. Hence, the coefficient of d^ In Equation (5.63) 
becomes: 

*P = i/P *» + * ‘V " ^ ~ b P ,5 7 
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However* the mass continuity restriction requires that: 


Cg_- ^ ^ Cs *= 0 


so that: 




(5,73) 


(5+74) 


The equations for the coefficients of Equation {5.63) are presented fn Table 
5-V, It should be noted that these coefficients are expressed as much as 
possible in a manner consistent with that used for the bounded skew-upwind 
differencing method, The mneumonic devices used In implementing QUD5 into 
TEACH differ From the current nomenclature because implementation of QUDS was 
initiated prior to the derivation of the equations for BSU05?, 

&,2,4 Boundary Conditions 


In guns, boundary conditions are applied directly to the finite-difference 
coefficients Tlstd in Table S-V. There are four types of boundary conditions 
used in the present analysis: 


[1) specified wall* 

[ 2 ] specified inlet* 

£3) a>tfs of symmetry 

(4) unspecified opening. 

The discussion presented in this section relies heavily oh remarks made 
earlier if) Section 5.1.4 for the bounded skew-upwind differencing procedure* 

5. 2.4.1 Scalars 

When calculating for the flu* on a face of a control voluine, adjustments are 
made if this face is near a boundary. For the west face, adjustments are made 
for the situations shown in Figure £-5. 





Case 1 : West face 
boundary 


Casa 2 : Wast-ses-t Face 
boundary, flow 
positive 


Casa 3 : East face boundary, 
flow negative 


Figure 5-9 Axial Velocity Boundary Condition Application 
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table i-v 


COEFFICIENTS OF THE FINITE-DIFFERENCE EQUATIONS FOR QUDS 


■s, - ^ (1 , ^ 

j V*w w V : v W 




XgT^W 


*E " + B E. + Ke-J!^ " fl ^= J 

- (1 - u > [] » ) -E^ 


X EE~\ 
^ET^K K EE“ 1( E 


Y -Y Y -Y 

P_ 5 ¥ N *SS 


“ ^ <W + % - C S , B U-o s ) ^ + "* c S « s IW v N .v s Y s -Vf 


- ”£ V„ <W 


V Y P V V P 
V Y S Y a- y s 


-‘Wri + % + ‘Vr, (1- B > 7^ - O-^) t*« n (1-0 ) Y ~ 7^ 

N s h s NN N 


Y -Y Y -Y 

n. Vv „ j- t , p s p s 

— fl — C„J t^Q,, (1”0 ! d) 

S ^ 5 S v -y V V 

M P CT S 


u , ''p-Jtw ^P'^HW 

*HK = -<V “w % 11 - “w> 

*B = ^ <1-“S S S 


V Y^-Ys Yp-Yg 

^ ‘ ” s “* U - “*> Ys-%S 


V i f j-Yp Y H -Yp, 

a - ft) On a - <*n) - — — - — ~ 

T ffiT‘P TflTT] 
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These boundaries can be any of the following typts; a) Specified opening, b) 
Unspecified opening, c) Specified wall, and d) Axis, Separate treatment is 
required for each of these types. In the fell owing sections, these boundaries 
are coral dared for each of the cases described above. 

Case 1 


If the boundary is a specified wall, axis, or an unspecified opening, the 
nodes W and MM are decoupled from the system of equations and, the influence of 
the will is incorporated into the finite- difference equations using the source 
terms; hence, Ay and Auy are set equal to zero. For a specified Inlet, the 
value of 4 is known and is effective at the west face, hence, Ay W is set 
equal to zero. However, because central differencing is used, Dy (the 
diffusion through the west facej. Is calculated as If jJ was located a full 
node away. Therefore, one extrapolates to get from fa t sip, and fa. 

The resultant extra terms are incorporated Into the source te^ms. 

Case 2 


The west face for the west node Is coincident with a boundary and the velocity 
on the west face Is eastward, in this case, the value of (Syy must be 
approximated. For an axis of symmetry, this value Should be the same as fa. 

For an unspecified opening, the assumption that the flow Is parallel to the 
gridlines In the region of the opening again suggests that P at the west- west 
node should be set equal to ^at the west node. For a specified wall, the 
node is set equal to Pu because for scalar variables the gradient 
In 3 near the wall Is assumed to he zero* For a specified opening, the known 
value on the boundary Is used, and an extrapolation can be used so that dyy 
Is a function of fay, and dp. In all cases, the resultant extra 
terms are incorporated Into the source terms. 

Case 3 

For flow westward through the west face when the east face is a boundary, the 
approximation for fa uses at the east node rather than on the east 
face. For an axis of symwetry, unspecified opening, or specified wall, the 
treatment is analogous to Case 2 and fa Is set equal to For a 
specified opening, an extrapolation of dy, fa* and tfe Is used to 
approximate fa. As In case 2* all resultant extra terms are incorporated 
Into the source terms. 

Since the computational molecule used in (JUDS is more extensive than that used 
for SUHS, there are a number of additional special cases that must be 
considered In modifying the coefficients of the finite- differ once equationsln 
the vicinity of the boundaries. These modifications are too tedious to examine 
here, and the interested reader should consult the computer program developed 
during this effort. 



5.3 COST OF QilDS AND ESUGS2 


In Table 5-VI „ the CPU time required to assemble the coefficients of the 
finite-difference equations relative to that required when using the hybrid 
differencing method are shown for B5UD5E and QUCS as implemented in the 
2D -TEACH computer code used in this contract. E5UDS2 requires considerably 
inore CPU time to assemble coefficients than does QUDS, since DUDS is 
unbounded* its coefficients can he calculated as needed fi.o*, "In line"). For 
BSUDS2, these coefficients must be computed and stored prior to application of 
the bounding strategy; then, the coefficients must be hounded in a process 
that requires nearly as much CPU time as that required for corputing then in 
the first Instance* Since the solver requires only about 30 percent of the 
total CPU tirse for any of the schemes, the total cost of using B5UD52 is 
nearly twice that of QUD5. These time requirements do not include the effect 
of the time required to achieve a converged solution using hybrid* bounded 
skew, or quadratic differencing. This effect varies from problem to problem. 

In some cases, 0UD5 required more iterations to achieve convergence; in other 
cases, 35UB52 required more iterations. Tn general d both E$UD$2 and i)UD$ 
required istore Iterations to converge than did hybrid differencing* In Section 
G detailed convergence Information for all test cases is given* 

The storage requirements for B^UDS-E and QUbS are also different {Table 5-fTI) r 
B SUDS 2 requires 20 extra two-dimensional storage arrays whore QUGS requires no 
additional arrays, in part because RSUDS2 was bounded whereas QUDS was not* 
However, 85UD52 cart also be programmed In such a way as to require less 
storage, but only at the cost of reduced speed and increased complexity of the 
code. Since 2D- TEACH Is a large code, the percentage increase In storage is 
net significant. 

In 30 -TEACH* BSUDS2 also requires tvfce as much CPU time as does hybrid 
differencing per Iteration* Since the number of coefficients which have to 
assembled In 3D are three tiroes as many as 2D-TEACH, It can be seen that 
3D- TEACH has been coded much more efficiently than has 2D -TEACH. The storage 
requirements of BSUDS2 are increased significantly when compared to 2P- TEACH; 
that is S3 additional 3D-storage arrays are required. Since 3D-TEACH analyzes 
less physically convex flows {apart from flow dimensionality) than 2D-TIACH 
analyzes, the Storage rsqui rement goes Up by a factor of two. 



TABLE 5- VI 


COST Of IMPROVED ACCURACY 


T T me 


Sc hems 

Coeff 1 cient Assembly Time 
(.1 teratTon x. node^ 

Solar 

[Iteration * node x swoop 

HYBRID 

l.o * (0.07m sec) 

1 + (0.015m sec) 

BUS052 

2.3 

1 

quds 

1*3 

1 


‘normalized to 

HYBRID 



TABLE 5 -VII 

COST OF IMPROVED ACCURACY 
Storage 

Humber of Estra 

Storage Increase 

Sc homo 

Arrays Used 

(percent) 

usuose 

m 

13 


0 

0 


$8 



6.0 DISCUSSION OF TEST CASES 

Two classes of test eases were used to demonstrate the accuracy improvement of 
the computer programs used fri this contract. The first class consists of a set 
of highly Idealized cases for which the exact results are known. This class 
was used 11) to verify that the computer programs were substantially free of 
programing errors and (£) to determine by numerical experiment the properties 
of selected mathematical features of the program. This first class of test 
cases was used In conjunction with the model problem studies and Is described 
In Section 4,0. The present discussion Is limited to describing the results 
obtained using the second class of test cases,. 

The need for the second class of test cases arises because In the first type 
of problems the flow field was assumed to be uniform and properties such as 
density and viscosity were constant. In addition, the mesh was uniform and 
cell aspect ratio was unity. Hone of these conditions exists when analyzing 
real combustor flows; hence, Conclusions reached after running the first class 
of problems can become invalid when real combustor flows are calculated, It Is 
therefore essential to Check the selected schemes In situations which ore 
simple enough in geometry to be modeled exactly by ED- TEACH and ap-TFACN and 1 
yet be representative of 3n -combustor calculations. In addition to Satisfying 
the above criterion, these test cases have to meet the conditions used to 
select bench mark experiments suitable for verifying CFB codes Iftef. 1 ) . These 
conditions are listed below, 

13 Minimum necessary flow dim ensionality. For the two dimensional version of 
TFAfH^ experiments in wh fcFT "the' "flows can be represented as two-dimensional 
were Selected. Three-dimensional flow situations were used only when 
Specifically testing the three- dittsensional version of the program. 

2 ) Veil -behaved fl ows. Flows in which instabilities! periodicity* or changes 
fn gross beha'vl o”f — Occur as flow conditions change were avoided* For example, 
flows were not used In which the location of reattachnent points of separated 
flow regions could undergo significant shifts as Soynolds number Is changed 
over the range of Interest. 

3) Known boundary conditions . Experiments were selected for which entrance 
[and exit) flow profiles could be specified as completely as possible. For 
turbulent flows* Initial profiles of turbulence intensity and integral length 
scale WEre sought. 

4) F xtensi Ve instrumentation . Attempts were nade to find flow mapping 
experiments in which noni’ntrusi ve techniques were used to characterize flows 
throughout the chamber volume. 

As the test cases selected originally were rgn using both the baseline and 
revised Versions of the computer program, it became apparent that additional 
test cases were needed. For example, the accurate prediction of turbulent 
flows requires the use of both a good numerical scheme and a sufficiently 
accurate turbulence model. Current turbulence models require further 
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development and their use may produce errors that cannot be separated from 
those produced by alternative numerical procedures. Therefore, two lanfnar 
flow cases were also selected for analysts; however, one of these cases was 
not based upon experimental data. 

The turbulent flow test cases were based on experiments discussed In 
References 1, 32 a and 33, In general, none Df the test cases satisfies 
completely ell of the criteria noted previously. In some cases, one or more 
Important Inlet profile Is not reported; In these cases, the unknown inlet 
conditions have been assumed or extrapolated from downstream results. However, 
all of the test cases do contain flow Features representative of those found 
In modern gas turbine engine combustors, 

Tho test cases used were: 

1) Laminar flow over a rearward facing step (E'R) 

£) Laminar flow with swirl In a sudden expansion [?D) 

3J Turhutent flow over a rearward facing step (2 d) 

4} Coannul ar turbulent flow in a sudden expansion [ ZD) 

5) toarwular turbulent swirling flow in a Sudden expansion (ED) 
fj) Row Of turbulent jets In turbulent crossflow (3D) 

6.1 TWO-DirHEMSlQHAL LAMINAR FLOW TEST CASES 

Two laminar flew cases were used. The first consisted of the flow downstream 
of a rearward facing step. The second consisted of the swirling flow 
downstream of a sudden expansion. These cases were run to assess the 
performance of the alternative finite-difference schemes independently of the 
turbulence nodal, 

fi.t.i Flow Downstream of Rearward Facing Step 

The laminar flow downstream of a rearward facing step was computed using the 
three finite-difference schemes considered herein; the baseline hybrid method 
- HYBRID, the second bounded skew- upwind differencing scheme ~ BS-UDS?, and the 
quadratic upwind differencing scheme - ODDS.. The flow geometry is shown in 
Figure 6-1, The Inlet Reynolds number was assumed to be Z5A, the Inlet axial 
velocity profile was taken as uniform^ and the inlet vertical velocity was 
assumed to be zero. The test case was run using four uniform, and successively 
finer, meshes designated TDPUAT1, TDPHATE, T0PHAT3 and TobKATfl. Mesh size 
parsrneters and computed flow reattachment points are presented in Table 6-1. 

It is seen that the mesh size is halved in each df recti on for each subsequent 
mesh. Representative computed streamlines are shown in Figure 6-Z. 
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TABLE 6-1 


CPMPUTFD ^ATTACHMENT LENGTHS FOR VARIOUS F1FFEREHCIHG SCHEMES 
FOR THE FIRST L Alii WAR FLOW TEST CASE 

Reattaehmeot Length, X R /h 


Case 

Grid 

HYBRID 

CUBS 

B5UDS2 

TOPIJATI 

10 x 6 

2,5 

5.32 

2.2 

T0PHAT2 

tD x 12 

3.14 

5.25 

3.35 

T0PHAT3 

40 x 24 

4.52 

5.35 

5.7B 

T0PHAT4 

78 x 48 

5,53 

5.75 

s.ge 


The variation of computed reattac foment length with total number of grid nodes 
{the product of the number of nodes In each direction} Is shown In Figure 6-3. 
Also shown Is the experimental ly determined value of 6.1 reported by Durst 
fRef. 3£) for this flow. It can be seen from Figure 6-3 that* for the coarsest 
nosh 5 HYP, RID and RSHDS2 have comparable accuracy that Is considerably worse 
than that of QUDS. The predicted reattachment length using QUCS Is not too 
sensitive to the number of grid nodes used. As the mesh Is refined, the 
accuracy of ASSJ0S2 increases rapidly and becomes somewhat superior to that of 
tjL'DS. It appears that all of the methods give essentially the same asyndetic 
value as ttie mesh is refined. For mesh densities normally used in 
two-diisensional flow calculations (i.e., "lODC nodes]. BSlTDS? and f)L'DS perform 
about equally weIT and both are superior to the HYBRID method. It should be 
noted that this assessment is specific to the type of flow calculated in this 
test case and is probably dependent on the Reynolds number of the flow. 
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NUMBER Of nodes; 

Figure o-3 Computed Reattachioent Length Using Different Finite Difference 
Schemes. 



From Figure 6-3, It can be seen that the computed reattachment point Is always 
less than the measured value, It fs possible that the underprediction of 
reattachnent length indicates that a grid- Independent result has not yet been 
achieved. More likely, this slight (~4 percent) error is due to the uniform 
axial velocity inlet profile used in the calculations since the experimental 
inlet profiles were not available from Reference 32. 

6*1.2 Swirling Flow Downstream of a Sudden Expansion 

The second laminar flow test case consisted of the swirling flow downstream of 
a sudden expansion in flow area. The Reynolds number for this flow was 450. 

The inlet tangential velocity was assumed to be equal to the axial velocity to 
give a vane angle of 45 degrees and a swirl number of 0*66* The expansion 
ratio, which Is defined as the ratio of the outer radius, R h to the inner 
radius, r, was 3.0. The flow *feld was terminated by a sudden contraction with 
a contraction ratio of 3.0. Uniform profiles of axial and tangential velocity 
were specified as the Inlet boundary conditions* There are no experimental 
measurements available for this case. The geometry and flow conditions are 
given in Figure 6-4. This test case was run using three uniform meshes 
designated COARSE, FINE, and XFINE. and described In Table 6-11, The computed 
streamlines shown in Figure 6-4 make it obvious that It Is difficult to 
characterize such flows using a single parameter such as reattachn*ent length. 
At a given station, the velocity profiles change significantly with a 
relatively insignificant change in the length of the recirculation zones. 



— 15r^ 

Figure 0-4 Flow In a Sudden Expansion with Swirl; Re & 450, Swirl Uo. => 
Lf.&6„ Vane Angle = 45-degree. 


TABLE 6- 1 1 

MESH DENSITIES FOR THE 
SECOND LAMINAR FLOW TEST EASE 

Case Srld 

COARSE 20 X 10 

FINE 40 x £0 


KFINE 


78 x 40 



Figure 6-5 shows the stagnation stream! fnes computed using the COAftSET mesh. It 
can be seen that the solutions produced by the three flnlte-df fferente schemes 
differ significantly from each other a although two basic reef relation zones 
are being calculated by all of the methods. The results for BSUDS2 show a 
distinct tobed central recirculation region while the results for GUOS- show 
that a lobe Is just beginning to form. The results for HYBRID shew no lobed 
region. For the FINE mesh, Figure G-6» hybrid differencing still does not show 
a lobed region hut detailed Inspection of the computed velocity profiles shows 
that further mesh refinement may lead to a lobed regfon. Computed stagnation 
streamlines using the finest mesh, XFlUE, are shown in figure 6-7, It can be 
seen that the shape of the recirculation region calculated using GSUDS2 has 
changed the least. The central recirculation region calculated using QUDS now 
shows a larger inner bubble which Indicates that for even finer meshes GUD5 
may also produce the- lobed recirculation bubble being computed by BSUtiSJ. 
■Results for H YE RI D still do not show break-up of the central recirculation 
region, but a comparison of the velocity profiles for FlWE and XFlKE meshes 
Indicates that further mesh refinement will lead to break-up. 
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Figure u-b Stagnation Streamlines for Second Laminar Flow Test Case on 
x 10 MeSti. 
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From the above, ft tan be concluded that the results far BSUDS2 are probe My 
nearer to a mesh- independent solution than the results using either HYSRJD or 
QLJDS, although ft appears that the use of QUdS will produce a mesh-independent 
solution more quickly than HYBRID. The implication of this conclusion is that, 
for this flow situation, the B5UD52 results appear to be more accurate than 
the QtfDS results. 

It Is Interesting to note that lobed recirculation iones have been observed in 
e>:per Indents for similar swirl numbers and geometries in turbul ent flow (ReF, 
31} as shown In Figure 6-8, It can be seen that ft is the""6acfi - pressure effect 
due to the downstream contraction that produces the Tobed feature* 
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Figure b-li Experimental ly Observed Lobed Recirculation Regions In Turbulent 
Flew - Yoon and Lilly (Ref* ] 


The results of the second laminar flow test case are In apparent contradiction, 
with those of the first test case, because in the first test case QUCS was 
■more accurate than RSUIJSJL Perhaps one explanation is the Haw angle L The flew 
angle in the first test case, Figure 6-2, seem* to be smaller than the flow 
angle of the second test case. Figure 6-4, It was found durfno the model 
Problem studies. Section 4.0, that the accuracy of these schemes was flow 
angle dependent and probably the flow angles of the above two cases are such 
that one favors PUDS and the other HSUDSi. 
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6 h £ TM -DIMENSIONAL TURBULENT FLOW TEST CASES 


Three two-dimensi anal turbulent flow test cases were run using the baseline 
and revised computer programs. The first test Case consisted of the analysis 
of the flow hewn stream of a rearward facing step. The second considered 1 the 
coannular nonswirling. axisyimnetrlc flow downstream of a sudden enlargement in 
flow area. The third case used swirling flow in the sane flow geometry as tihe 
second case. 

fi.£.l Turbulent Flow Downstream of a Rearward Facing Step 

The first test case consisted of the analysis of the flow downstream of a 
rearward -fad rig step. Fk tensive data have been reported for this configuration 
In inference 33. This experiment has been used by several workers to assess 
the accuracy of their flow field models and it has been designated as Stanford 
Case m. 

The geocsetry of the test section used in the experiment reported in Reference 
32 Is shown in Figure 6-?; the vertical scale has been exaggerated for 
clarity. The flow region from a point ^ step heights upstream of the step 
(where suitable data for initial conditions were obtained} to a point 20 step 
heights downstream of the step (well beyond the location of the reattachment 
point) was analyzed. In Figure 6-D t this region is delineated by the dashed 
lines, 
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Figure Geometry of Test Section for Flow Over a Rearward-Facing Step - 

First Turbulent Flow Test Case 


Initial conditions for the baseline case are shown in Table 6- HI and are 
based on test conditions reported in Reference j 4. The fritfal boundary layer 
was approximated In a step-wise manner as indicated in Table S-III. 
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TABLE G-II1 


INITIAL CONDITIONS F0fi CASE m 


Mean Velocity 13-32 m/sec 

Temperature 41.2 C 

Pressure 1 01 - 35 KM/m*- 

boundary Layer Thickness 9-4 ram 

Turbulence Intensity Eu'/u) .W3 


Axial Velocity Profile - Assumed 


u * 9.16 d/sk, 
y ■ 13.74 m/sec, 
u = 18.32 m/sec * 


1.0 s y/h <1.0336 
1.0336 sy/h <1.0736 
1.0736 <y/h ^3.0 


The turbulence energy dissipation rate was calculated from the expression 

t m 3.0K 3/Z /H 

where K is the turbulence kinetic energy and V is the height of the channel. 
The expert men tally -determined reattac Indent point was located at 7 + 1 
step-heights fron the step. 

Three cases, each using a different grid system, were calculated using the 
baseline computer program for the conditions shown in Table 6-1 n. The 
pertinent grid information Is given In Table 6-V and the grid set up for three 
cases is explained below, The first case used the relatively coarse grid shown 
in Figure 6-1(3 and the grid has been designated as COARSE; in Figure 6-10, the 
vertical scale has been exaggerated for clarity. A second case was then run 
using a grid obtained by nearly doubling the number of grid lines in each 
direction; this grid has been designated as FINE. Cased on the results 
Obtained for these two cases, a third case was run In which the number of 
vertical grid lines was approximately the same as that used In FIVE, but the 
distribution of both horizontal and vertical grid lines was changed to reduce 
the Peel et number distribution ip the vicinity of both the step and the 
reattacliraent point; this grid has been designated as ADJUSTED. In the baseline 
code, the finite-di f fere rice approximation to the governing flow equations uses 
central differencing for Peclet numbers less than 2 and upwind differencing 
for Peclet numbers greater than Z. Errors due to numerical diffusion are 
associated with the use of upwind differencing. The Peclet number distribution, 
for the horizontal (east-west) direction for ADJUSTED Is shown In Figure 6-11. 
It can be seen that, even for this dense mesh, the Peclet number Is in the 
range from 2 to ID for a relatively large portion of the separated flow 
region. Since the streamlines are at an angle to the grid In this region* 
numerical diffusion is Still present. fThe high Peclet number zone, PeeTet > 
10, outside of this region has little influence on the computed results 
because the streamlines are almost parallel to the mesh, ) The Peclet number 
distribution for the vertical (north-south] direction is shown in Figure 6-12 
and indicates that central differencing Is used for the convective terms in 
the vertical direction throughout the separation zone* thus, the calculated 
results should be relatively free of errors due to the finite-difference 
approximation for these terms. 
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Figure u-H* Llrid System for Stanford Case 42H -COAftSE* 


The tabulated and measured reattachment points are presented in Table 6-IY. 
Use of either of the alternative grid systems (that is t fine or ADJUSTFD) 
resulted In a substantial improvement in the accuracy of the predicted 
location, As the results in Table 6-IY shew, ADJUSTED is reasonably free from 
numerical diffusion. Further mesh refinement is required to make the 
calculation completely free of numerfeaT diffusion. but the number of nodes 
needed to eliminate diffusion was beyond the storage capability of the 
computer used in this program. 


TABLE 6-1V 

CALCULATED AND MEASURED ftEATTACHHFNT POINTS 
fh = step height) 


Measured (mean value) Xp/h ■ 7 


Cal tula ted: 

Case 

HYBRID 

Roattachmem length „ 
BSUTI52 

. V* 

jjjJDS 


COARSE 

5.2 

£.4 

5*5 


FINE 

5*7 

5.9 

Unstable 


ADJUSTED 

5,8 

5.8 

Unstable 
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Flyure &-11 East’Hest Peclet Nurtier Distribution for Stanford Case 

421-AajjtfSTELi Using HYURip Differencing {Pe * Peclet Humber) 



AXIAL DISTANCE, */h 


Note; [X and Y axis are not on same scale) 

Figure e-12 tJortr^Sduthi Pectet Number Distribution for Stanford Case 

421 -ADJUSTED Using HYBRID Differencing (Pe ■ Peclet Humber} 
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F or completeness, additional information about the grid system used arid the 
number of Iterations required to achieve convergence is summarized in Table 
6-Y* 

TABLE 6-Y 

ADOlTttm GRID SYSTEM INFORMATION 



Humber of Grid Lines 

Number of Iterations to Convergence 

Case 

K-D1 recti on Y-Direction 

RiB Rip 

BUSES? 

COARSE 

26 

29 

303 

253 

FINE 

5fl 

56 

7BD 

533 

ADJUSTED 

74 

S3 

1515 

993 


The first test case was then run using the BSUDS2 finite-difference method , 

The grids and flow conditions were identical to those used with the baseline 
computer program using hybrid differencing. In fact, the same Input files were 
used to run both the HYBRID and B5L1D52 cases. The calculated reattachment 
points are presented In Table 6-IV where it can be seen that BSJD52 is 
5 lightly more accurate than hybrid differencing for both the COARSE and FINE 
meshes. However * the improvement In accuracy Is not as large as suggested by 
the model problem studies* Section 4.0* and laminar test case 1. For the 
ADJUSTED mesh* the reatfachment point given by HYBRID and BS0DS2 Is the sane. 
It is curious to note that the reattachment length with the ADJUSTED grid 
decreases slightly as compared to the FINE grid. 

Since the AWUSTTD grid was designed for HYBRID to give more accurate results* 
it is not surprising that this grid has become somewhat Inefficient for 
fiSUDSE. The number of iterations to achieve convergence is less for the BSDDS2 
than for the HYBRID method. However* since PSUDS2 requires the calculation of 
more finite- difference coefficients than HYBRID, total computation time nay be 
greater when DSHD 52 is used. 

Results of running the above test case with DUD'S are also given fn Table 6-1 V. 
For the coarse grid* it can be seen that the trend shown in laminar test case 

1 is followed. QUDS is marginally better than RSUPS? and HYBRID. However, the 

difference fn accuracy is not a; large as it was in the laminar flow test 
cases. For finer grids, DUDS becones unstable. This behavior was expected, 
because as dfscussed in Section 4.0* £)UDS is not compatible with the solver 
presently incorporated in TEACH. 

ft. 2,2 Eoannular Monswlrllng Turbulent Flow 

The computer programs were applied to the analysis of the coaxial nonsulrling 
Flow downstream of a sudden enlargement in flow area, as reported by Johnson 
and hennett in Reference 35. The flow is admitted into a 51 nm radius tube by 

two concentric annuli. The inner annulus has an exit radius of 1£.S rm„ and 

the outer annulus has an exit radius of approximately 15.3 mm. The 
•experimental apparatus Is shown in Figure 6-13, The flow in both annuli Is 
water. 
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Figure 6-13 Second Test Case - Co-axial jets with artd without swirl 


The tube dividing the two annuli is tapered to a sharp edge at the exit of the 
annular region (that Is* at the entrance of the sudden enlargement region ) . 
Very limited data arc currently available at the end of the annular region* 
hence, the foil owing assumptions were made to establish the boundary 
conditions for the Inlet* 

Up stream of the annular jet* measured profiles of axial velocity were 
provided. Turbulence energy was calculated from measured turbulence intensity 
assuming Isotropic turbulence. Since no measurements were available for the 
length scale, the following expression was used for the calculation of 
dissipation rate' 


c . ( 3,0 K^l/Oh 

where % Is the hydraulic diameter* 

£x peridental Information on the boundary conditions for the Central jet wore 
provided at a station 12,7 mm downstream of the expansion. Measurements could 
not toe made at the expansion plane Itself due to lack of optical access* It 
was therefore decided to use the Information provided at 12*7 mm to estimate 
the conditions at the sudden expansion* in the following manner: 

o Axial Velocity - From the measurement of peak velocity at the 12*7 mm 
station and the flow tfeynolds number, a Fully developed turbulent 
pipe flow profile was calculated which was then adjusted to give the 
correct mass flow rate. 
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o Turbulence Energy - The measured value of the turbulence energy on 
the centerline at 15.7 nr, normalized by the centerline axfaT 
velocity at this point, was used to provide a Uniting boundary value 
of turbulence energy at the inlet* assuming that this ratio remained 
constant: 


K/u 2 • D.005 

Turbulence energy dissipation rate was calculated fn the same manner as ft was 
calculated for the annular jet* 

Global values of the Initial conditions are presented in Table 6-lfI* The axlaT 
velocity profile for the inner annulus was specified in accordance with the 
velocity distribution calculated using the analyses mentioned above* 

TABLE G-VI 


INITIAL CONDITIONS FOR NOKSWIRLIWG FLOW TEST CASE 
(Based on Reference 35? 


Mean Velocity* nv/sec 
Mass. Flow Bate, fcg/sec 
Temperature, “C 
Turbulence Intensity (u'/u) 


Inne r Flow 

0.596 
0* 64fl 
15.6 
0,057 


Outer Flew 

1,74 

0*55 

15*5 

0.036 


Three cases, each using a different grid system, were calculated for the 
conditions shown In Table 6-VI* The iirst case used the relatively fine grid 
shown In Figure 6-14 and has been designated EXTRA FINE grid. The second case 
used a grid with the same vertical definition as that shown in Figure 6-14 but 
with approximately one-half the number of axial nodes; this case has been 
designated FINE grid. The third case used approximately pre-half the number of 
nodes in each direction* relative to the grid shown id Figure 6-14 and has 
been designated COARSE grid. The vertical scale in Figure 6-14 has been 
exaggerated for clarity* 

A compari son of the measured and calculated centerline axial velocity 
distribution is shown In Figure S-1S for the coarse grid. It can be seen that 
in the initial region the results do not agree well with the data regardless 
of the differencing scheme used. 

Figure 6-16 shows the comparison for the FINE grid. Since &HJS was unstable 
for this mesh, only HYBRID and £$UD££ calculations are shown. It can be seen 
that there is virtually no difference between the two calculations. However, 
it cannot be said that 9 mesh Independent solution has been obtained, because 
east-west cell FficTet numbers are large even for this grid* Figure 6-17* 

Worth- south cell Pec let numbers. Figure 6- IB, are well within limits* 
Calculations with the EXTRA FINE grid did not show any further Improvement and 
hence are not shown* 
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Figure u-16 Axial Velocity uHtrfbutloji - Second Turbulent Flow Test C$se 



Figure 6-17 East-West Peclet Humber Distribution fur Second Turbulent Flow 
Test Case; Pe ■» Peclet Number. 
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Figure 0-1 tl North- South' Peclet Number Distribution for Second Turbulent Flow 
Test Case; Pe c Peclet Number* 


Figure 6 - 1 & shows the calculated streamlines for the case* This figure also 
shows the three stations at which comparisons of the radial profiles of the 
mixture fraction and axial velocity were made. These comparisons are shown fn 
Figures 6-20 and 6-21 for the coarse grid and in Figures and fi-23 for the 
FINE grid. It can be seen that the agreement with the experiment is reasonable 
even for the coarse yrid, and nestr refinement does not result in further 
improvements. The apparent good agreement is somewhat misleading because 
comparison of Che centerline profiles had shown quite a significant 
discrepancy between the data and calculations. Also* as was found in the first 
turbulent flow test case* Che accuracy of hybrid differencing fs almost as 
good as that of BSUDS2 and QUDS* and emphasizes the need to investigate 
further this behavior. 
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The calculated and measured reattachment points are shown: in- Table 6-Y]|. The 
measured reattachment point was estimated from flow visualization data. It can 
be seen that the trend fs the same as was seen for laminar flew over a 
rearward facing step* Figure 6-3. B5UDS2 seems to reach an asymptotic value 
more quickly than does hybrid differencing. But, as pointed out earlier the 
difference between BSUDS2 and hybrid differencing is not large. Since QUfiS 
became unstable* its accuracy cannot he judged for this test case. 


TABLE 6-Y11 

COMPARISON OF REATTACHMENT LENGTH FOR THE CORNER 
RECIRCULATION ZONES FOR COANNULAR NQNSHlRLtNG FLOW 
MEASURED 200-22 5mm 


KE$H CALCULATED RE ATTACHMENT LENGTH 


HYBRID 

DSHD S 2 

ms 

mm 

mm 

mm 


COARSE 

2 JR 

241 

234 

FINE 

243 

254 

UNSTABLE 

EXTRA FINE 

252 

257 

UNSTABLE 


For completeness i additional information about the grid system used and the 
number of iterations required to achieve convergence is surmari led in Table 

g-yiii* 

TABLE 6- VI 1 1 

ADDITIONAL GRID ST STEM INFORMATION 


Case 

Grid 

Iterations to Convergence' 

k 

HYBRID 


pirns 

COARSE 

21 x 20 

367 

300 

1911 

FINE 

21 x 38 

13.95 

1228 

Unstable 

EXTRA FINE 

40 x 38 

15D6 

1430 

Unstable 

* Convergence defined as 

maximum residual 

source (normalized) 

0 r DQE 


6*2.3 Coannul ar Swirling Turbulent Flow 

The test case for swirling flow is based upon the data reported by Johnson and 
Etobacfc {Ref h *nj ) for the experimental apparatus which was used for the 
coannul ar nonswirling flow* An extensive number of comparisons between 
experimental and analytical results have been reported previously under the 
HDST/Aerothermal Modelling Program - Phase I (Ref. I) using the hybrid 
differencing scheme. 
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Since very limited data are currently available at the end of the annular 
region, the fell owing assumptions were made to establish the boundary 
conditions for the Inlet In the swirling flow case. The calculations were 
started at the first plane of measurement* 0,51 cm from the expansion plane* 
using the measured quantities as Inlet boundary conditions, At this plane all 
variables* except the dissipation rate* required for the 2D-TEACH code were 
measured and can bE specified, To specify the dissipation rate, it was assumed 
that the inner portion of the flow originated from the inner tube and the 
outer portion originated from the concentric outer tube; turbulence length 
scales were assigned accordingly* using 3 percent of the passage heights, ami 
the dissipation rate was calculated from these length scales. 

Four cases* each using a different grid system, were calculated using the same 
set of inlet conditions. The first case used the relatively coarse grid shown 
in Figure and has been designated COARSE grfd test case. The second case, 
designated ADJUSTED COARSE grid, used the same number of grid lines as coarse 
grid but with a greater concentration of grid lines in the axial direction 
near the entrance region. The third case* designated FIIJE grid, used a greater 
number of grid lines in the axial direction than either coarse grid or 
adjusted coarse grid; again* the grid lines were concentrated near tho 
entrance region. The fourth case used the same axial grid Tine distribution as 
FIKE grid together with more grid lines in the radial direction; this case has 
been designated EXTRA FINE grid. EXTRA FlWE grid is shown in Figure 5-25, The 
number of grid lines for all casts in both the axial and radial directions, 
together with the number of iterations necessary to achieve a converged 
Solution, are listed In Table G-IX, 
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Figure d-24 Grid System for Coarse Grid Test Case 
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Figure i)~£s> briti System for Extra Fine Grid Tost Case 


TABLE 6- IX 

COMNLLAR SWIRLING TURBULENT FLOW 
ADDITIONAL GRID SYSTEM AND CONVERGENCE INFORMATION 





Iterations to Convergence 


Case 

Grid 

HYBRID 

BSUDS2 

DUOS 

COARSE 

<0 x 35 

955 

11G6 

3616 

COARSE ADJUSTED 

40 x 35 

2594 

1340 

1579 

FINE 

47 x 35 

1539 

1422 

1240 

EXTRA fine 

47 x 43 

2221 

364* 

1676 


*Run with FAST Algorithm, which is the PM version of PISO f Pressure-ImpTfcIt 
Split Operator) Algorithm - Kef. 40 , 
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Contour plots for the east -west Peclet number di strfbuticns obtained using the 
baseline -computer program are shown for COARSE grid and EXTRA FINE grid cases 
In Figures 6-26 and 6-Z7, respectively. There were very few differences 
between the ADJUSTED COARSE, FINE and EXTRA FINE results; hence, only the 
EXTRA FINE grid results are being shown. These contours show the region for 
which the absolute value of the Peclet number Is ID or less. Genera- ly* the 
hybrid differencing procedure can be used for Peclet numbers less than ID, but 
In flow regions with strong gradients in flow properties such as the entrance 
region and corner recirculation region In the present case, the Peclet numbers 
should be less than 2 to ensure the use of central differencing and resulting 
solution accuracy. Inspection of Figures 6-26 and 6-Z? shows that these 
criteria arc not met for even the finest grid. However, It can be seen that, 
as the grid is refined in the a*1al direction, the east-west Peclet numbers 
are reduced to more acceptable values In both the entrance region and the 
corner rec i nc ul atl on region. 
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Figure u-2ti East-West Peclet Number Distribution for Coannular Swirling Flow 
- Coarse Grid 
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Figure o-27 Easi-West Peclet Number Distribution for Cearmular Swirling Flow 
- Extra Fine lin'd 


Contour plots for the north-south Peclet number distributions are shown in 
Figures 6-28 and 6-29, Generally, the radial velocity components are small 
relative to the axial velocity components, and, therefore, the north-south 
Peclet numbers are smaller than the east-west values. It is interesting to 
note that north -south Peclet numbers change significantly from the cbarsE grid 
to EXTRA FIME grid. Figure 6-2 8 and 6-29, The reason for the change in 
north- south Peclet numbers is that increasing the grid density in the axial 
direction modifies the flow field significantly, producing much largEr radial 
velocities. These increased radial velocities are reflected in the change in 
north- south Peclet numbers. The north- south Peclet numbers, as was the case 
for east-west Peclet numbers, do not change significantly for the last three 
grids* and their magnitude is such that the grid In the radial direction is 
probably adequate. 

Before detailed profile comparisons are made for this flow, it is instructive 
to compare thE predicted reattachment lengths of the corner and central 
recirculation regions. These lengths are defined as the distance between the 
Inlet and reattachment point Figure 6-30, The length of the corner 
recirculation region is given in Table 6-X. It can be seen that for hybrid 
differencing and B5UDS2 the predicted length is reduced by almost a factor of 
two between the coarse and adjusted coarse grids; further refinement does nut 
bring about a significant change In. this length. However, for [HUDS, it can be 
seen that this length is reduced only by about ID percent between the coarse 
and EXTRA FINE grid. Since QUflS seems to reach an asymptotic value sooner than 
does B5UDS2 and hybrid differencing, the implication is that for this case 
QUDS seems to be more accurate than B$DD$2 and hybrid differencing. 
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Fiyure u-30 streaml i r*e$ for Coanrtul ar Swirling Flow 


TABLE 6-X 

COMPARISON OF REATTACHfCHT LENGTH FOR THE CORKER 
ftEC I RCULATIQN REGION FOR COAMNEJLAR NON SWELLING FLOW 


MESH SCHEME 


hrBRID 

B5UD5? 

ms 

mm 

\m 

rnri 


COARSE 

113 

ID 2 

sa 

COARSE ADJUSTED 

54 

56 

50 

FIRE 

49 

54 

51 

EJiTRA FTNE 

50 

53 

52 


Table 6-XI gives the reattachment length for the central reci irul ation regions, 
and again {RJD5 seems to be attaining an asymptotic value more quickly than 
BSUDS2 or hybrid differencing. In fact, hybrid differencing does not approach 
an asymptotic value even for the finest mesh* The superiority of PUDS is 
confirmed when the centerline velocity profiles , Figures. 6-31 and 6-32, and 
radial profiles of axial and tangential velocities, Figures 6-33 to 6-36, are 
examined. It can be seen that for the coarse grid, CRJDS agrees best with the 
experimental results* When tho mesh is refined, there is no significant 
improvement in QUDS; however, both HYBRID and B5UD-S2 short significant 
Improvement and come close to QUDS results. It was also found that 11511052 was 
unstable for the finest mesh; art Improved pressure algorithm, FAST, was used 
to achieve convergence* It Is possible that QUDS has reached grid Independence 
because Its results did not change when the mesh was refined. However, the 
grid refinement was not sufficient to state this with confidence. 




TABLE 6-XI 


COMPARISOH OF REATTACHWENT LENGTH OF THE CENTRAL 
RECIRCULATION REGION FOR COAHNULAR NONSUIftLING FLOW 


MESH 

HVBRID 


mT' 

COARSE 

120 

COARSE ADJUSTED 

135 

FINE 

130 

EXTRA FINE 

142 


SCHEME 


IISUDS2 

QUSS 

im 

mm 

134 

136 

146 

155 

153 

15S 

156 

157 



Figure o-31 Axial Velocity Distribution - Third Turbulent Flew Test Ca 
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Ffyure o-32 A*i al Velocity distribution - Third Turbulent Flow Test Case 
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Figure L-J3 Comparison of Axial Velocity Profiles - Coarse Grid 
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Figure b-34 Comparison of' T&nijential Velocity Profiles - Coarse Grid 


SWI1LEHG CDAHtolLAR TURR-ilLENT TLDV 

HYBRID 

GUD5 - — — 

BSUDS2 



Figure u-3i Comparison of Axial Velocity Profiles - Extra Fine Grid 
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Figure o-Ju Uomparison of Tangential: Velocity Profiles - Extra Fine Grid 


6*3 SUMMARY QF 2D -TEST CASES 

Two laminar flow test cases and three turbulent flow teat cases were 
calculated anri. compared against experimental data where such data were 
available* For laminar flow, the results confirmed the findings of the model 
problem studies; both schemes were markedly superior to the hybrid 
differencing scheme, BSUDSZ seemed to he more accurate at large flow angles, 
whereas QUCS was more accurate at small angles. For turbulent flows, 
especially for the first two test cases, hybrid differencing turned out to be 
almost as accurate as the two improved accuracy schemes. This behavior was 
unexpected and warrants further investigation. Also* QUDS became unstable for 
the first and the second turbulent flow test cases. This instability of tJLrhS. 
was expected* For the third turbulent flow test case* QUPS was superior to 
both B5LfB$2 and hybrid differencing for the coarse mesh. It was found that 
BSUDS2 was unstable for the third turbulent flow test case and the improved 
pressure algorithm, FAST, was used to achieve convergence. This behavior of 
b$l)t)5£ was also unexpected and needs further Investigation. The unexpected 
behavior of the differencing schemes in turbulent flow emphasizes the 
Importance of using realistic test cases for the final evaluation of the 
differencing schemes* 

6.4 SELECTION QF DIFFERENCING GCHEHE FOR 3D-TEACH 

It was found from the computation of ZB- test cases that was prone to 
Instability. Since in 3D* this tendency is bound to increase, there was little 
choice but to select the more stable B SUDS 2 for incorporation Into 3D -TEACH* 

6.5 THREE-DIMENSIONAL TURBULENT FLOW TEST CASE 

The baseline f hybrid differencing! and BSHDS2 versions of the three- 
dimensional computer program were used to analyze the three-dimensional test 
cases. This test case is based on the experiments conducted by Khan {Ref. 37! 
in which rows of jets were injected normal to the main flow in a duct with a 
rectangular trots- section. For the test condition simulated* a single row of 
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£.54- cm diameter jets with a pitch (jet center! Ins- tc-centerllne distance) to 
jet diameter rat To of four was used; the test section height to jet diameter 
ratio also was four. A schematic diagram of the flow geometry In the test 
section Is shown in Figure 6-37, 


z 



Figure u-37 Test (iase Flow Geometry for Base-Line Computer Program. 


In the axial tx) direction, the simulation extends from a plane five jet 
diameters upstream of the leading edge Of the jet to- 0 point 24 jet diameters 
downstream of the trailing edge of the jet; that Is, f&r 30 jet diameters. In 
the vertical direction, the calculation extends from the floor of the wind! 
tunnel to the roof, a distance of four jet diameters* In the lateral direction 
the computational domain covers the distance through the jet centerline to the 
plane of symmetry between the Jets, a distance of two jet diameters* Kahn 
verified tl)at the flow was symmetric about the latter plane. 

The grid system for eaOh plane Is shown In Figures 6-3S through 6-40, 
respectively* The jot was simulated using the relatively fine, rectangular 
grid in the Injection plane, as shown schematically in Figure 6-41, The number 
of grid nodes in the x, y, and t directions were 34, IQ, and 16, respectively. 
For later reference, this grid system has been designated as the COARSE grid. 
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Figure o-Jtf Hoarse Mesh for the Test Case fn X-Y Plane 
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Figure &-41 Grid definition in Vicinity of Jet Exit. 


The ratio of jot to mainstream velocity was approximately 2,3, Experimental 
data, Reference 38, were used to establish the initial velocity of the jet. 
These data, which were for a jet to mainstream velocity ratio of ?, are shown 
in Figure 6-4£, together with the assumed jet velocity di stribotfon used in 
the test case. For other calculationsfn this flow system, see Ref. I 

For the baseline computer program with hybrid differencing, calculations were 
continued: for 141 Iterations until the imaxim-un" residual error was less than 5 
percent * Selected calculations were continued to a lower level of maximum 
residual ; the computed flow field did not change significantly. 
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Figure b-42 Jet Axial Velocity Profile, 


Both BSUDS2 and Skew-Upwind Differencing ( SUD > without bounding wore used* The 
number of iterations required to achieve comparable levels of convergence were 
as follows: 


Differencing 

Hybrid 

BSUDS2 

SUD (Ho Bounding) 


No* of Iteration 

141 

111 

244 


The B51D52 calculations were run with the blending factor Initialized to unity. 


Calculation time per iteration for BSUD5Z was approximately twice that for 
hybrid differencing. This large Increase Is probably due to the time required 
to assemble the additional finite- difference coefficients required for BSUH52 
and to the bounding procedure, as noted earlier* 

The streak line plots for the coarse mesh are shown fn Figure 6-43 to 6-45, 
These figures do not show significant differences between the hybrid and 
BSUDS2 calculations* 
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Figure i>-44 Comparison of Streakin' ne Plots for B5UD52 and HYBRID in the Y-Z 
Plane at Y/D ■ O.ij - Coarse Grid 
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Figure Comparison of Streakline Plots for B5UD52 afld HYBRID in the X-Z 

Plano at Y/EJ = 20 - Course Grid 
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Comparisons of calculated results vizi) data at four a x f a 1 locations are shown 
In Figures 6-46 to 6-49. As expected* calculations without bounding appear to 
contain less numerical diffusion than the bounded calculations, However,, the 
difference between the two predictions is not significant. Also* either skew 
method Is successful; in reducing numerical diffusion* as can he seen from 
Figure 6-46 where peaked profiles are calculated by the skew schemes where as 
numerical; diffusion smears out the peaks for hybrid. There are* however, 
significant differences between the experimental data and the predictions, and 
grid refinement is required to ascertain If the discrepancy Is due to the 
coarseness of tha mesh or Inaccuracies of the turbulence model. 
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Figure u-46 Comparison of Calculated and Measured Axial Velocity Profiles at 
X/b “ 4 from Jet Centerline and 2/0*0 Using Coarse Grid 
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Figure ti-44 Comparison of Calculated and Measured Axfal Velocity Profiles at 
K/D * B from Jet Centerline and Z /D = 0 Using Coarse Grid 
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Figure Comparison of Calculated and Measured Axial Velocity Profiles at 

= 1 ti froffl Jet Centerline and 1 /D = 0 Using Coarse Grid 



The mesh used far the fine grid calculations used 40, 20, and 17 nodes fa the 
*, y z direction* respectively, (Figures fi“5Q to fi-55) it can be seen that the 
major difference between the coarse and fine meshes is the grid density fa the 
vertical direction, wihlch has doubled. 
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For FliJE mesh calculations, the I>1 ending factor was initialized to 0.7&. The 
reason for initializing the blending factor to 0.75 was that, when this factor 
was initialized to 1.0, (as is usually the caseh severe convergence problems 
were encountered. These problems persisted even when a converged hybrid 
differencing solution was used as the initial guess for the BSLfDSZ s cheap. 
Since the initialized blending factor is not altered at a given node if the 
solution remains within bounds, reducing the initial blending factor hy 2S 
percent nas the effect of introducing some extra numerical diffusion. This 
diffusion seems to be sufficient to damp out the oscillations and to enable 
tne solution- to converge. Increasing the blending factor beyond 0.75 made the 
Solution unstable. It was later found that the computer program was not 
computing the blending factors properly in some cases and the above 
instability could be due to the error fn the code. Since it was ensured that 
the blending factors were In bound, the final solution Is not in error. Jt may 
have more than the minimum amount of numerical diffusion required by the 
BSUH32 scheint!, However, as will be seen later, it has significantly less 
numerical diffusion than hjurid. 

yQkflparisujis of calculated results with data are shown in the Figures 6-5 1 to 
It can be seen that as compared to the CU/HSE grid calculation, the 
predicted profile shapes have been modified, end the hybrid differencing and 
BSUU52 predictions approach each others however, there are still significant 
differences between the predictions and the data. Examination of the 
streamline plots* Figures 6-57 to 6-59, confirms the findings of the profile 
plots. The differences between the flow fields predicted by the two schemes 
are less pronounced than they were for the coarse grid. 
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Figure 6-53 Comparison of Streamline Plots for B SU D S 2 and HYBRID in the X-Z 
Plane at ijfy =■ O.S - Tine Grid 
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Figure d-e9 Comparison of Streamline Piets for BSHDS? and HY&FtlD fn the Y-Z 
Plane at X/D - £0 - Fine Grid 
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Figure d-bO shows a comparison of COARSE and FINE mesh predictions for B5UDS2 
and hybrid differencing., It can be seen that the predicted profiles have 
changed siyni flcantly near the floor of the tunnel but not near the roof. A 
comparison of the CGAftSE and F1KE meshes in the y-z plane* Figure 6-61, shows 
that although the mesh was refined considerably near the floor its density was 
not increased sufficiently near the roof. Since the measured profile near the 
roof is also peaked and the gradients near the roof are as Targe as those near 
the fleer* it -can be Inferred that a significant amount of numerical diffusion 
Is still present in the FINE grid calculation. To ensure that a mesh 
Independent solution has been achieved* the mesh density near the roof must be 
increased^ Unfortunately* the cost of making such a run during o resent 
contract was prohibitive. 
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Figure o-6{> 3D Test Case 1 - Comparison of Anial Velocity Profiles for -Coarse 
and Fine lies he s 
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Ft guns ts-— 6 1 JD Teat Case 1 - Campari son of Coarse and Fine Meshes in the f-Z 
Pi ane 


tt.G SLfr^lAftV OF 3D-TEST CASE 

Due turbulent flow three-dimensional test case* a row of Jets in cross How, 
was calculated using hybrid differencing and 0SUU52. Two mesh sites, COARSE 
P4 x 10 x 15 ) and FINE (40 x £0 jq 17} were used. It was found that BSUDS2 was 
able to reduce numerical diffusion when compared to hybrid differencing. There 
remained a considerable discrepancy between the data ami predictions. Because 
of budget limitations, a finer mesh, which would have teen able to resolve the 
velocity near the roof of the tunnel, could not he used. Hence* it could not 
be ascertained if this discrepancy is due to the turbulence inode 1 or to an 
insufficient number of grid nodes. 

0.7 SUMMARY OF TE5T CASES 

The improved accuracy f ini te-di f fe rente schemes, QUDS and &SUDS2* were 
evaluated In 2D-TEACH by computing laminar and turbulent flow test cases, 

Tnese computations were compared with experimental data where such data were 
available. Every test case was calculated using successively finer grids in an 
effort to obtain grid Independent solutions. Grid Independence was obtained in 
most cases* but in some turbulent flow test cases* limitations of computer 
size or instability in the solution prevented achievement of this goal. 




For the laminar Flow test cases, the flow ever a backward facing step and 
swirling flow into a sudden expansion were analyzed. For laminar flow, the 
results confirmed the predictions of the model problem studies described in 
Section 4,0. Both schemes were markedly superior to the hybrid, differencing 
scheme; anti, BSUDS? Was more accurate at large flow angles and QdOS was more 
accurate at small angles. Three turbulent flow test cases were considered: 
flow over a backward-facing step, nonswirling flow in a sudden expansion, and 
swirling flow in a sudden expansion. For the first two test cases the accuracy 
of hybrid differencing was the same as that of the two improved accuracy 
senemes. This behavior was unexpected and warrants further investigation, 

Also, QUDS became unstable for the first and second turbulent flow test cases. 
This instability of QUDS was expected, BSUDS2 became unstable for the third 
turbulent How test case; this behavior was also unexpected and needs further 
investigation. The unexpected behavior of the performance of the differencing 
schemes for the turbulent flow test cases demonstrated the inadequacy of 
rely in g only upon model problem studies for selecting a scheme and emphasizes 
tne Importance of using realistic test cases. 

One three dimensional* turbulent flow test case, a row of jets in cross flow, 
was calculated using hybrid differencing and B5UU5? with two mesh sites. 
Because of computer storage and budget 1 imitations, a grid-independent 
solution could not be achieved for either scheme. It was found that E5UD52 was 
more accurate than hybrid differencing for this test case* 
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7,0 CONCLUDING REMARKS 


TmIs stud/ has contributed significantly to the effort of reducing run&rfcal 
urrer in computational fluid mechanics cedes. One of the most 'Important 
(level opuunts has been the establishment of a frame work by which future 
schutaes can ce judged. This fratse work did not exist prior to this study, In 
this stuuy, not only were tlie criteria for evaluation established, hut a 
nui.ioer of scheites were evaluated using these criteria and the procedure was 
successful In Identifying the most suitable scheme for 3D-eodes. In addition, 
tiid assessment wtJiods developed during this study were able to predict most 
of tne deficiencies of tno selected scheme which were liter exhibited In the 
Ju-code, Another fmportant outcome of this study was the realization that an 
improved accuracy scherae cannot be completely divorced from the solution 
alguritiM* It was because of the Incompatibility of QUDS with the present 
solution algorithm. and the resulting instability, that it was not 
incorporated Into db-TEACn* although it has significant advantages over BSUD52. 

Tins effort nas resulted In the Improvement of the accuracy of 2G as well as 
jb-cqdes baseo on the THAtK concept. Although accuracy improvement is 
dependent on the flow field, ft Is expected that for practical applications, 
unere s-wne portion of the flow makes a large angle with the grid Tines, this 
improvement will be more than sufficient to justify tne increased cost. Since 
t»e selected scneme is never less accurate than hybrid differencing, there Is 
nv danger of the improved code giving less accurate answers than the baseline 
version, 

Tnis study should be looked upon as i first step in the process of reducing 
numerical error from combustor performance codes. More work is required, 
especially for the 3b code, to ensure that the selected scheme is performing 
with, optimum accuracy, 

Tne reasons for tno poor performance of USUDS2 in swie of the 2D- test cases 
also need to be Investigated more throughly. At present, it Is conjectured 
t.iat trie sensitivity of each scheme to flow angle is responsible for this 
ueliavior. It is gdite possible that this may not be tike only reason. For 
example, all model problem studies were conducted or meshes of uniform size 
with cell aspect ratio of unity whereas the 2b turbulent flow test cases used 
a non-uni form mesh uitn large aspect ratio cells, hence poor performance of 
8 £Ul£ 2 in ill case 5 could also be attributed to large cell aspect ratio. In 
principle, modifications to the scheme can be made to take account of the cell 
aspect ratio, The scheme can also be modified to account for the pressure 
gradient iri the flow Field which, as discussed earlier, leads to inaccuracies. 
It is suggested that tine test cases be rerun with the above modifications to 
tne scheme to assess tne Improvement in performance. It is anticipated that 
the performance of H5UDS2 will improve significantly after these modifications. 

Tne ul end lug procedure used in BSUDS2 suffers from the disadvantage that the 
final solution Is hot path independent. It is even conceivable that different 
initial values of the blending factor may produce different Convergedsolutions 
wiiere convergence is uaasured by the level of residual source error which in 
practice is always nonzero; that is, the solution nay not be unique 
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CorjputatiDiidl 1/, Path dependency may cause difficulties when this code is used 
In a production cnvIroKiflieiit, ievurul sc^&me-s can be formulated to reduce or 
eliminate this path dependence. One scheme would be to always start t^lth a 
converged .Hybrid solution, The other would L>e to recalculate t^e l>1 ending' 
Factors at every Iteration and use their latest value without under relaxing 
trie™ with tueir value of the previous iteration. 



e.L> RECOMMENDATIONS 


There are two type of recommendations that can be made based on this work, 
tfecoimneudaticns of the first type are related to the use of current codes for 
optimum accuracy, kecojimendatlofis of the second kind are concerned with 
furtner developments of the code. 

It is recommended that &SUUS2 be used for production running of £0- and 
3P-TEACK codes. Since flows conmonly encountered: to gas turbine combustors are 
complex* and since some portion of the flow [makes a large angle with the grid, 
lines, E5U0$£ will yield more accurate results than will hybrid differencing. 
It Is also recoin, lended that during grid selection account be taken of property 
gradients present in the flow field* Indt sertmt nate refinement is not feasible 
in three-dimensional calculations with present day computers. 

Since the accuracy QUUS as well as that of &SUB5S is dependent on flow angle* 
It is desirable to develop differencing schemes that are less sensitive to the 
direction of the velocity vector. It is not necessary to abandon present 
methods completely in the pursuit of flow angle independence. Since QUDS and 
SUES are accurate at different flow angles. It may be feasible to combine 
tnese two schemes in such a manner that the resulting scheme vs more accurate 
than either of the schemes at all flow angles. This possibility should be 
investigated. It was found that the model problem based an single cell 
calculations was Successful In predicting the flow angle sensitivity of the 
Scheioes tested. Hence, it is recommended that this problem be used to screen 
tiie candidate schemes in future, 

Incorporation of tne more implicit pressure algorithm* FAST, into the 3D-TEAEK 
code is recommended * It Is expected that the stability of BSUli$2 will Improve 
with this algorithm as demonst rated in the 2D-tode. In order to improve the 
stability of QUDS and other schemes, whose computational molecule is not 
compatible with the AUI method* a more compatible matrix solver must be 
developed. 

Alternatively, instead of developing further the 3B-TEADH code, an entirely 
new computer program based (say) on a tl-me-marching procedure which does not 
need a matrix solver. Since the time -step size necessary for stability of an 
explicit algorithm is extremely small , it is expected that a higher-order 
di fferendng scheme, such as quas, will be stable In this framework* However* 
Considering that the time required to bring a new code to production Status Is 
about ten years and the expenditure of I'esobrces is considerable, any 
improvements to the present code In the interim are highly desirable* 
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APPENDIX A - THE THREE DIMENSIONAL BOUNDED SKEW-UPWIND DIFFERENCING SCHEME 


In this section, the three dimensional Bounded Skew-Upwind Differencing Scheme 
(BSUDS) is described. First, a brief review of the flux form of the equations of mo- 
tion is presented. Second, a detailed description of the finite-difference form of 
the flux contribution to a representative face of a typical control volume is given; 
the derivation of the flux contributions to the other five faces is then outlined. 

Third, the resulting co-efficients for the finite-difference equations representing 
the total flux (and sources) are presented. Fourth, the results of applying the boun- 
dary conditions in a manner consistent with the foregoing flux representation are 
shown. Fifth, the bounding scheme for the coefficients is detailed. 

This section is an extension of the two-dimensional scheme presented in Section 5.0 
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- Flux Form of the Equations of Motion 


The equations of motion for both laminar flow and (time-averaged) turbulent flow 
can be written in similar fashion for all of the dependent variables: 

Cartesian Co-ordinates 



+ 1 3 

r\ ' r\ 


- — (I i\ i!) + s 

r 96 r ‘P 96 <1 


The variable d represents any of the dependent variables (e.g., the velocity compo- 
nents u, v, w, mixture fraction, turbulent kinetic energy and turbulent energy dissi- 
pation). The exchange coefficient, T^, represents the sum of both laminar and turbu- 
lent contributions and is interpreted as the effective viscosity for (j> = u, v, w, the 
effective diffusivity for cf> = mixture fraction, etc. is a generalized source term. 

Eqs . (Al.l) are integrated over a control volume appropriate for each dependent 
variable 4> and, after some manipulation, the finite-difference equivalent forms of 
Eqs. (Al.l) are obtained. The control volumes are defined using an orthogonal grid 
formed by the intersection of co-ordinate lines in each axial, radial, and lateral 
co-ordinate direction. The intersection of the grid lines form the grid nodes at 
which all flow properties except the axial (u) , radial (v), and lateral (w) veloci- 
ties are calculated; i.e., all scalars. The axial velocity is calculated using a 
second grid with grid nodes located midway between the scalar grid nodes in the axial 
direction and co-incident with the scalar grid nodes in both radial and lateral posi- 
tions. The radial velocity is calculated using a third grid with grid nodes located 
midway between the scalar grid nodes in the radial direction and co-incident with the 
scalar grid nodes in both axial and lateral position. The lateral velocity (or azi- 
muthal velocity in cylindrical co-ordinates) is calculated using a fourth grid with 
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grid nodes located midway between the scalar grid nodes in the lateral (z or 6) direc- 
tion and co-incident with the scalar grid nodes in both axial and radial position. 
Directions in the grids are identified as north, south, east, west, back and front. 
Nodes are identified as shown in Fig. A-l where P is the node at the center of the 
control volume. For hybrid (upwind and central) differencing, the principal nodes are 
designated N, S, E, W, B, F and P. For the bounded skew-upwind differencing scheme, 
twenty additional nodes are used as shown in Fig. A-l. 

The relative positions of the various grid systems can be illustrated by con- 
sidering the grid systems in the x-y plane only. Then, the grid system for the scalars 
is shown in Fig. A-2 and the control volume for the scalars is shown in Fig. A-3 . The 
relative positions of the axial and radial velocity control volumes is depicted in 
Fig. A-4. Similarly, in the y-z plane, the radial and lateral velocity control volumes 
are displaced relative to the scalar control volumes. Finally, in the x-z plane, the 
axial and lateral velocity control volumes are displaced relative to the scalar control 
volumes . 

The faces of the control volumes for each scalar are defined by planes located 
midway between the scalar grid nodes as shown in Fig. A-3. Thus, the location of the 
east or west faces of a scalar control volume is co-incident with the axial locations 
of the axial velocity components, etc., so that normal velocity components lie along 
the boundaries of the control volumes. (Generally, boundaries for the u, v or w con- 
trol volumes include scalar grid-lines and are not necessarily located midway between 
velocity grid lines.) For the scalars, the grid systems defined above provide some 
computational convenience. Since the u, v and w velocities are stored midway between 
the scalar grid nodes, the convective fluxes for each face can be calculated without 
recourse to averaging any of these velocities. Also, the pressure gradient driving 
the flow can be computed without averaging pressures. 

The finite-difference forms of Eqs. (Al.l) are derived by integrating these equa- 
tions over the appropriate control volume. In performing the integration over the con- 
trol volume for each term in Eq. (Al.l), the mean-value theorem is employed and the 
source term is linearized in the vicinity of the center of the control volume (point 
P). After some manipulation, the finite-difference form of Eqs. (Al.l) is obtained 
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CE^e 


where C , 
h 


where a , 
e 

"diffusion 


C W^ w + C N^n “ c s * s + C B <t> b " C F*f 

= D„ (d> -d> ) - D r , (<b -d> ) 

E v P W vv P 

+ \ W - D s W 

+ D b <VV - D F W * s 0 + Vp 


: TT , etc. are "convective coefficients" as defined below 

w 

C! = (pu) a 
E e e 

C W = (pu) w a w 

C N = (pv) n a n 

C S = (p v ^ s a s 

C B = (pw) b % 

C F = (Pw) f a f 


(A1.2) 


(A1.3) 


i . a , a , a, , a _ are the areas of the faces of the control volumes. The 
w n s b f 

coef f icients"are given by 




(A1.4) 
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where the convective and diffusive coefficients are shown in terms of the cartesian 
co-ordinate system. Since Eq. (A1.2) has the same form for both cartesian and cylin- 
drical co-ordinate^ it is only necessary to define the geometric parameters used in 
the convective and diffusive fluxes in a manner appropriate for each coordinate 

system. These definitions are given in Appendix B. If a consistent set of 

geometric parameters is used, then the equations for the bounded skew-upwind diffe- 
rencing scheme developed below are identical in both cartesian and cylindrical co- v 
ordinates. For convenience, the nomenclature for cartesian co-ordinates will be used. 

It is important to note that Eq. (A1.2) applies to all of the dependent variables 
although the appropriate grid must be used in each case to define the geometric 
parameters used in the calculation. Also, Eq. (A1.2) can be used with any of the difference 
procedures considered in the program since each scheme is simply an alternative method 
for defining (interpolating for) the fluxes at the face of the control volume (e.g., 

<j) , $ , <j> , cf> , Au 4> )• However, the diffusion terms are always represented by 
e w n s TO f 

central differences. 

It is convenient to define a total flux for each face of the control volume as 
the sum of a convective flux and a diffusive flux such that 

F e " F w + F n - F s + F b ~ F f = S u + (A1.5) 

where 

F e = C E^e " D E (<t>E " <t>p) 

F w = Cvj<f> w - % ( p " <t>w) 

F n = C N^n % (< *>N " 

F s = c s^s “ D S 

F b = C B^b “ d B ~ ^ 

F f = C F^f “ D f (<}> P " M (A1.6) 

In the following section, the skew differencing procedure will be used to calculate 
the values of the dependent variables at the faces of the control volume. As a result, 

Eq. (A1.5) will include not only the values of <J> at the "normal," or main, grid node 
locations (E, W, N, S> B, F and P) but also at the corner locations. The finite-difference 
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form for Eq. (A1.5) is then 


Vp = Ve + Vw + Vn + Vs + Vb + Vf + ^NE 

+ A SE^SE + ^W^NW + A SW*SW * A BNE^BNE + "^BSE^BSE 
+ ^NW^BNW + a bsw^'bsw + A BN^BN + A BS^BS + a be^be 
+ A BW^BW + a fne^fne + A FSE^FSE + A FNW <J) FNW + ^SW^FSW 
+ A F N*FN + WfS + A FE^FE + A F W^FW + S u + Vp 


(A1.7) 
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A. 2 - Calculation of the Fluxes 


Recall that the equations of motion, Eqs. (Al.l), can be written in terms of fluxes 

to each face of the control volume, Eq. (A1.5). In this section, a procedure will be 

described to calculate fluxes, F , F , F , Fu and F.e. The derivation of F„, the flux 

e w n D r w 

to the west face, of a typical scalar control volume is given in detail. The deriva- 
tion of the other fluxes and of the fluxes of the u, v, and w velocity components are 
outlined . 

Consider the control volume used to determine the value of a scalar variable at 
the point P. A portion of this control volume is shown in Fig. A-5 for the case in 
which the u and v components of velocity are non-negative and the w component is nega- 
tive. For all velocity vectors located at the center of the west face of the control 
volume (point w) , the flux to the west face is given by: 

F «"V,- D » ( *p-V <A2 - 1) 

Furthermore, for central-differencing (CD), the value of the dependent variable at the 

west face, <j> , is always given by linear interpolation between $ at the W and P grid 
w 

nodes : 


w 


= (i - 


a ) 

w 


W 


+ a 


w 


where the interpolating factor is 


(A2.2) 


a - X w ~ X W 
w X P “ X W 


(A2.3) 


In the 3D-TEACH computer program, a = 0.5 for each face of the control volume for 

each scalar variable since the control volume faces are located midway between scalar grid 

nodes. For the axial, radial and lateral velocity components, however, the a for 

each face may assume other values. The central difference form of the flux at the 

west face is then: 




a ) + 1 
w 




(“w^e^ 


■ ~ 1 ) 
w P 


where the Peclet number at this face is given by: 


(A2.4) 
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(A2.5) 


_ w 


= C/D. 


W W 


The central difference form of the flux at all faces of the control volume for point P 
are presented in Table A-l . 

The upwind difference (UD) form for the flux at the west face is defined by: 



5 w W 


(A2.6) 


for a non-negative value of axial component of velocity, u^, and 


F 

W UD 

D W 


V <f> 

e w v P 


(A2.7) 


for a negative value of u^. Essentially, the upwind differences are derived from the 
central difference by setting a w equal to zero or unity, respectively, and by neglect- 
ing diffusion. 

It is convenient to define the velocity switches 


° u -M 1 + 17r) 

„v. % ( 1+ ^) 


for the axial, radial, and lateral velocity components and to adopt the convention that 
the value of a switch is unity if the velocity component is non-negative and it is zero 
otherwise. Then, the upwind difference form of the flux at the west face of the 
control voluem is: 


Fw UD 

= P e 

D W 


o u (f>„ + (1 - o u ) <J>. 
w w w 1 


(A2.9) 


The upwind difference forms at all of the faces are summarized in Table A-2. 

Equation (A2.1) is also the starting point for developing the flux equations for 
the skew-upwind differencing (SUD) procedure. The value of <j> at the west face of 


167 



the control volume, <j> , is determined by extrapolating the velocity vector upstream 
of the face. For the case shown in Fig. A-5, this point lies in the plane formed by 
the nodes W, SW, BW and BSW. The value of w'" is determined by linear interpolation. 
First, the projection of the velocity vector in the x-y plane is extrapolated upstream 
to the point w* which lies along the grid line between the nodes W and SW. The inter- 
polation factor k-j is defined by: 


k 


1 


(w 1 W) = 1 Ax ^w 

Ay 2 Ay U w 


(A2.10) 


where (w'W) is the distance between point w' and the node W and Ay is the distance 
between nodes W and SW. Then, 


^w' k l (<i! SW ^ + ^W 


= k l <S> sw + (i - k-| ) <p 


1 J 9 W 


(A2.ll) 


For very large flow angles (skewing) relative to the co-ordinate directions, k^ will 
exceed unity and w' will be defined in terms of <j> at the SW and S nodes; however; it 
is known that for two-dimensional flows this approach can yield negative coefficients 
at corner nodes which can in turn produce oscillations in the solution; it is assumed 
that the same situation will prevail for three-dimensional flows. To assure that the 
coefficients for the corner nodes are non-negative, then: 


k^ = min 




(A2.12) 


The use of absolute values in Eq. (A2.12) permits this equation to be used to define 
k^ for all velocity components at the west face. Similarly, the value of cj> at the 
point w" along the gridline connecting nodes BW and BSW is given by: 

4>w" = k l ^BSW + (1 “ k x ) <t> BW (A2.13) 


where the same interpolation factor is used as in Eq. (A2.ll). The value of <j> at the 
point w'" is obtained by linear interpolation between 4> w ' and The interpolation 

factor is obtained by extrapolating the projection of the velocity vector in the x-z 
plane to the point 3 along the gridline connecting nodes W and BW: 
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(A2.14) 


++— 

* w - = k 


3 ^w" 


(1 - k 3 ) 


f w 


with 


kg *= min I l s % 


Ax 
A z 


w 


w 


The superscript ++- indicates U w :> 0, v w > 0, W w < 0. 
Consequently, 


V = k 3 [ k l ^BSW + ^ " k l> $BW ] 


(A2.15) 


+ (1 - kg) [ kg <j> sw + (1 - kg) <j> w ] (A2.16) 

For two-dimensional flow, kg = 0 and it can be seen that Eq. (A2.16) has the correct 
limiting value. 

Equation (A2.16) was derived for the case of U w >_ 0, V w 0, W w < 0. For the 
case of U w > 0, V w > 0 and W w >_ 0, the appropriate equation can be derived by rotating 
the velocity vector into the positive z direction and (referring to Fig. A-l) replacing 
the subscripts BSW by FSW and BW by FW. As a result, one obtains: 

= kg [ kg <J> FSW + (1 ~ kg) $ FW ] 

+ (T - kg) [ kg 4>g W + (1 - kg) <j>y ] (A2.17) 

Then, for U w > 0, V w > 0, and all W w 

= c” 4> wllI [ Eq. (A2.17) ] + (1 - a") [ Eq. (A2.16) ] (A2.18) 

or 

= o” kg [ kg <j> FSW + (1 - k^) <b F W 1 

+ (1 - a W ) kg E kg (JjbSW + ^ ” kg) 4>BW ] 
w 

+ (1 - kg) [ kg <j>SW + ^ - kg) <J>^ ] (A2.19) 


where the superscript ++ indicate that U w and V w are both non-negative. 


169 



A similar expression can be derived from Eq. (A2.19) for the case U w 0, 

V w < 0, and all W w by replacing the subscripts FSW by FNW, BSW by BNW, and SW by NW 
to obtain 

k 3 ^ k l ^FNW + (1 " k l 5 *FW ] 

+ (1 c”) k 3 [ k x <j> BNW + (1 - k x ) 0 BW ] 


+ (1 - k 3 ) [ k i ^ + (1 - k x ) ^ ] 


(A2.20) 


Then, one can combine Eqs . (A2.19) and (A2.20) to obtain an expression for cp ^ for 

U > 0, all V , all W from: 
u w — w’ w 


,+ _ v ,++ 

t* ,ii' 

w w w 


+ (1 - o ) 
w 




(A2.21) 


For the case in which U„ < 0, V„ > 0, and all W„, an expression for <t> can be 
obtained from Eq. A(2.19) by replacing the subscripts as follows 


Hence 


FSW becomes FS 
FW becomes F 
W becomes P 
SW becomes S 
BW becomes B 
BSW becomes BS 


4>~m " o* k 3 [ k i <f> FS + (1 - k l> <f> F ] 

+ (1 - o W ) k 3 [ k x <t> BS + (1 - kj) <p B ] 


+ (1 - k 3 ) C kj_ <j> + (1 - ki> p p ] (A2.22) 

Finally, for the case of u w < °> v w < 0 and W w» an expression for 4> can be 
obtained from Eq. (A2.20) by replacing the subscripts as follows: 

FNW becomes FN 
FW becomes F 
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BNW becomes BN 



BW becomes B 


NW becomes N 
W becomes P 


Hence 

V.. = °W k 3 [ k l ^FN + (1 " k l } *F ] 

+ (1 - 0») k 3 [ k x 4 > bn + (1 - k x ) 4 B ] 

+ (1 - k 3 ) [ kj_ <j) N + (1 - k l ) <j, p ] 

Then for U w <_ 0, all V w and all W w 

_ v v, 

4,m ° 4, + ' 0 ' 4.. in 

WWW w 


(A2.23) 


(A2.24) 


and for all U w , all V w and all W w : 


V" 


u 

a 

w 



+ (1 - a U ) 
w 



(A2.25) 


In analogous fashion, interpolated values for 4 at the other faces of the control 
volume can be derived. Equations for <p ,,, , <p ,,, , and 4^,,, can be derived immediately 
from the expressions for 4 w m , 4 tTI » and 4^ m , respectively, by shifting nodal sub- 
scripts in each case in the positive co-ordinate direction (i.e., toward the opposite 
face of the control volume). Tables A-4 through A- 9 present the equations correspond- 
ing to Eqs . (A2.19) through (A2.25) for each of the faces. 

Since the grid is not uniform, the geometric parameters ax, A A z > nsed in the 
calculation of k^ and k 3 differ for each face of the control volume and for each 
orientation of the velocity vector at each face. A summary of the geometric parameters 
is presented in Table A-3. Only the definitions for the parameters used for the east, 
north and back faces are shown in Table A-3. Since the fluxes at only these faces 
are actually computed as indicated in the next section; the flux at the west face for 
the control volume at node (I, J, K) is the flux at the east face of node (1-1, J, K) , 
etc . 
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In his original development of the skew upwind differencing procedure for two- 
dimensional flow, Raithby(Ref» A-1 )assumed that the value of the dependent variable 
at (say) the west face of the control volume is equal to the interpolated value. In 
the present development, this is equivalent to assuming that <t> = 4> ,,, so that — using 

Eq. (A1.6) and the definition of Peclet number— 

^WoTTr) 

-5— ■ p e w V" ‘ “W ( *P - *W> (A2 ' 26) 

w 

It is desirable to use the central-difference procedure for small values of the grid 
Peclet number and the skew upwind differencing method for large values of the grid 
Peclet number. It is also desirable that these two formulations produce a continuous 
transition at the transition Peclet number which in the case of U w > 0 at the west 
face of the control volume becomes: 


P* = — (A2.27) 

w a w 

For the scalar grid system used in the 2D-TEACH computer program, P g * =2. At the 
transition Peclet number, the central difference result using (Eq. A2.4) is: 


Fw CP _ 

I>W a w 


(A2.28) 


while the skew upwind differencing method (Eq. (A2.26)) yields: 


Fw SUP = K'" 
D W a w 



(A2.29) 


From Eqs . (A2.18) through (A2.25), it is clear that these two results are not equal. 

✓ ++ \ 

(Note: At present, 4> w m = <J> . ) 

The fluxes at the transition Peclet number can be made equal by noting (contrary 
to the assumption made by Raithby) that <j> w and 4>^,„ are related by: 


<}> 

Y w 


in 




As + . . . 


(A2.30) 


so that Eq. (A2.26) becomes 
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Fw SUP 

% 


= p 





(d> -(h) 

VH P V 


(A2.31) 


where As is defined in Fig. A-5. Writing the central differencing result in terms 
of the flux definition, Eq. A2.1, then: 


Fw CP _ 

Dy " Pe w 4> w (<i) P V 


(A2.32) 


Clearly, these two fluxes will be equal at the transition Peclet number if a correc- 
"k ( 3(j) \ 

tion Pgl J As is added to the skew upwind differencing flux, Eq. (A2.26), to 
obtain: 


■ w 


e,.,'*',.,'" T t e w ( If ) As - (<!)? - 4> w ) 
w 


SUD ^ * 

V = P e.A,., - P 


(A2.33) 


The derivative ( 3<|>/ 3 s ) w can be computed by 

ab ^w _ V" 


3s 


As 


(A2.34) 


To provide a continuous variation in flux at the transition Peclet number, the derivative 
in Eq. (A2.34) is approximated further as: 


3j; 

3s 


(<f>w) CD _ ♦w'" 
As 


(A2.35) 


where ( <f>w^ CD ^- s gi ven b y Eq. (A2.2). Then, after using Eqs. (A2.2) and (A2.35) and some 
rearrangement, Eq. (A2.23) becomes: 


F 


w 


SUD 


W 


(Pe w ~ P^) 4> w «" + C Pe" (1 - « w ) + 1 ] 4 W 

+ («*w p e w " 1) <f> P 


(A2.36) 


This last result is valid for all U , V w , W w (i.e., for all <j> w »n) provided that the 
transition Peclet number is given by: 


p* = 1 

e w o u — ( 1 - a ) 
w w 


(A2.37) 
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At the transition Peclet number, the flux given by both central differencing and skew 
upwind differencing are: 



Results similar to Eqs. (A2.36) and (A2.37) for all faces of the control volume are 
summarized in Table A-10. 

In the hybrid differencing procedure, the more accurate central differencing formu- 
lation (see Table A-l) is used when the Peclet number is less than the transition value 
while the less accurate, but stable, upwind differencing result (see Table A-2) is 
used when the Peclet number is greater than the transition value. In the development 
of the equations for two-dimensional flow, it was believed at first that a similar 
hybrid differencing method could be formulated using central differencing for P e < P e 
and skew-upwind differencing for P e > P e " . However, this approach proved to be 
unworkable since some of the coefficients derived from this hybrid formulation could 
be negative. As an alternative, a flux blending scheme was used in which a weighted 
average of the upwind differencing fluxes is used. The weighting (blending) factor, 

Y, was chosen to assure boundedness (i.e., all co-efficients of the difference equa- 
tions are non-negative). It is assumed that a similar flux blending procedure can 
be utilized for the three-dimensional case. The Bounded Skew-Upwind Differencing 
(BSUD) scheme is defined by 


Y n F n.c 


+ <1 ~ Y n> F i 


(A2.40) 


where II represents one of the six faces of the control volume. The weighting factor 
Yjj is restricted to the range 0 <_ Yjj ^ 1 and is calculated as described in Section A. 5. 
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A. 3 - Calculation of the Co-efficients for the 


Finite-Difference Form of the Equations of Motion 


The finite-difference form of the equations of motion (e.g., Eq. (A1.7)), can be 
derived directly from the flux information presented in Tables A-l , A-2, A-4 through 
A-10 and Eq. (A2.40). The resulting expressions will contain the unknown blending 
factor, y . The blending strategy requires that the terms in the equations for the 

coefficients responsible for producing negative co-efficients be isolated so that 
appropriate values for y can be determined. Furthermore, the coefficients for the 

control volumes adjacent to the physical boundaries of the flow may require modifica- 
tion to incorporate the effect of the boundary conditions. Thus, to simplify manipu- 
lation and modification, some additional notation will be defined. 


Let the center of the control volume (point P) be located at the Ith axial posi- 
tion, Jth radial position, and Kth lateral position. The flux contributions (the com- 
ponents of the total flux) to the east face are denoted as El(I,J,K), E2(I,J,K), 
E3(l,J,K,L); the flux contributions to the north face are denoted as Nl(l,J,K), 
N2(I,J,K), N3(I,J,K,L), and the flux contributions at the back face are denoted as 
B1(I,J,K), B2(I,J,K), B3(I,J,K,L) where L = 1,2, 3, 4 as indicated below. The flux con- 
tributions are defined as follows: 

_Centra_l Dif f erencing 


ei(i 

s J ,K) = 

d e “ 

ct C-p 
e E 

E2(I 

,J,K) = 

El(l. 

,J,K) 

E3(I 

,J,K,L) 

= 0 

L 

N1(I 

,J,K) = 

% " 

“n C N 

N2(I 

,J,K) = 

N1 (I ; 

,J,K) 

N3(I 

,J,K,L) 

= 0 

L 

B1(I 

,J,K) = 

% - 

°b C B 

B2(l 

,J,K) = 

Bl (I 

,J,K) 

B3(I 

, J,K,L) 

= 0 

L 


C E 

1 , 2 , 3, 4 


1 , 2 , 3, 4 


1 , 2 , 3, 4 


(A3.1) 


(A3. 2) 


(A3. 3) 
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Bounded jakew-Upwind Differencing 


E1(I,J,K-) = (o“ - 1) C E 
E2(I,J,K) = E1(I,J,K) + C £ 

E3(l,J,K,l) = D E (Pe e -Pe*)(k 3 k 1 ) e 

E3(I,J,K,2) = D E (Pe e -Pe*) [ k 3 (1 - k x ) ] g 

E3(I,J,K,3) = D E (Pe e -P e ;> [ (1 - k 3 ) k x ] e 

E3(I,J,K,4) = D E (Pe e -Pe*) [ (1 - k 3 ) (1 - k x ) - 1 ] g (A3. 4) 

Nl(l,J,K) = (oJJ - 1) C N 

N2(I, J,K) = Nl(l ,J ,K) + C N 

N3(I, JjK,l) = D N (Pe n -Pe n ) Oc 3 k 1 ) n 

N3(l, J,K,2) = D N (Pe n -Pe^) [ k 3 (1 - k x ) ] ^ 

N3(I,J,K,3) = D N (Pe n -Pe*) [ (1 - k 3 ) kj_ ] n 

N3(I, J,K,4) = %(?e n -Pe^) [ (1 - k 3 ) <1 - kj_) - 1 ] ^ (A3. 5) 

bi(i,j,k) = (o* - i) c B 

B2(l, J,K) = B1(I,J,K) + C B 

•k 

B3(l , J ,K, 1) = D B (Pe b -Pe b ) (k 3 k 1 ) b 

☆ 

B3(I,J,K,2) = D B (Pe b -Pe b ) [ k 3 (1 - k x > ] b 
B3(l s J,K,3) = D B (Pe b -Pe b ) [ (1 - k 3 > kj_ ] b 

B3(I , J,K,4) = D fi (Pe b -Pe b ) [ (1 - k-j) (1 - k x ) - 1 ] (A3. 6) 

The use of central vs. bounded skew-upwind differencing is determined by the value 
of the Peclet number at each face. The parameters, C^,, a e> k b , . . . are local 
values; the subscripts (l,J,K) have been omitted in the interest of readability. 
The corresponding flux contributions at the west face are given immediately by 
El(l-1 , J,K) , E2(I-1,J,K), ...» the flux contributions at the south face are 
Nl(l,J-l,K), N2(I,J-1,K) .... and the flux contributions at the front face are 
B1(I,J,K-1), B2(I,J,K-1) , . . . 
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The coefficients of the finite-difference form of the equations of motion may 
then be defined in terms of these flux contributions. The results are presented in 
Table A-ll. The notation used in Table A-ll can be explained by examining, for ex- 
ample, the second term, exclusive of the blending factor, in the equation for A g at 

the node (l,J,K): [ c u (1 - o w ) E3 ] (I,J,K,2). The subscripts I,J,K indicate the 
u w 

node at which a , o and E3 are evaluated; the fourth subscript refers to the flux 
e e 

contribution E3(I,J,K,L) with L equal to 2 — see Eq. (A3.1) or Eq. (A3. 4). 

For the two-dimensional case, it is relatively easy to demonstrate that Ap is 
equal to the sum of the other eight co-efficients when use is made of the mass con- 
tinuity equation. It is assumed in the three-dimensional case that A is now equal 
to the sum of the other 26 co-efficients in Table A-ll when use is made of the mass 
continuity equation: 


C E - 


c w + c 


N 


c s + 


C B ' 


C F = 0 


(A3. 7) 


The boundary conditions (Section A. 4) and the blending scheme (Section A. 5) can be 
applied directly to the flux contributions so that the results shown in Table A-ll 
are general. 
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A. 4 - Boundary Conditions 


The boundary conditions are applied to the flux contributions as defined by Eqs . 
(A3.1) through (A3. 6) so that the set of coefficients (e.g., Table A-ll) for each de- 
pendent variable is the same throughout the computational domain. This consistency 
simplifies the application of both (1) the algorithm for solving the set of simultane- 
ous equations for each variable and (2) the bounding procedure (see Section A. 5). 

_Scalars 

It will be recalled that the axial velocities are stored midway between the 
scalar gridlines in the axial direction; therefore, the axial velocities are stored at 
the midpoint of the east and west faces of the scalar control volumes. Similarly 
the radial velocities are stored at the midpoint of the north and south faces of the 
scalar control volumes and the lateral velocities are stored at the midpoint of the 
front and back faces of the scalar control volumes. 

Physical boundaries are located midway between scalar gridlines in the appropriate 
direction. In the following discussion, the boundary conditions are applied to the 
various faces of the control volume for an interior node (I,J,K) to produce modified 
east-west flux contributions (El , E2, E3) , north-south flux contributions (Nl, N2, 

N3), or back-front flux contributions (Bl, B2, B3) for use in calculating the co- 
efficients of the finite-difference equations for the scalars (e.g., Table A-ll). 

There are five types of boundary conditions to be considered; these types are: 

(1) specified wall, 

(2) axis of symmetry, planes of symmetry, 

(3) unspecified opening, 

(4) specified inlet, and 

(5) periodicity. 

The application of these boundary conditions to modify the flux contributions used to 
determine the finite-difference coefficients for a node (l,J,K) having a physical 
boundary on the west face of its control volume will be described in detail. The 
modifications to the co-efficients due to the other physical boundaries are similar 
and will be summarized. 
a) West Face 

For a specified wall, the value of the scalar at the wall may be either known (e.g. 
turbulence kinetic energy = 0) or unknown (e.g., turbulence dissipation). Since the 



variation of the scalar from its wall value to the value at the node (l,J) is 
(generally) non-linear, it is prudent in every case to decouple the wall values from 
the system of equations. The influence of the wall can be reflected by using appro- 
priate wall functions to compute the source terms for the scalars. From Table A-ll, 
it can be seen that the (known or unknown) wall values of the scalar can be excluded 


if the co-efficients A ,A , A , A ,A ,A,A, 

BNW BSW BW’ FNW FSW’ FW’ NW’ 


A and A are 
W SW 


These coeffiicents will vanish if the following flux contributions are set 


set to zero, 
accordingly : 


E2(I-1,J,K) = 0 

E3(I-1,J,K,L) =0 L = 1,2, 3,4 

N3( I , J ,K,L) = 6 
N3(I , J-l ,K,L) = 0 

’ L = 1 2 

B3(I,J,K,L) = 0 

B3(I,J,K-1,L) =0 (A4.1) 


If the west face of the control volume for the node (I,J,K) is a plane of symmetry, 
then the flux normal to this face is zero. By definition, the unknown values of the 
scalar at the nodes (l-l,J,K) and (l,J,K) are equal; the gradient of the scalar normal 
to the axis vanishes there. Thus, an axis of symmetry is mathematically identical to 
a specified wall boundary condition with zero gradient (source) at the wall. The modi- 
fications to the flux contributions are given by Eq. (A4.1). 

For an unspecified opening, it is assumed that (1) the flow is exiting through 
the opening, (2) the streamlines are parallel in the vicinity of the opening, and (3) 
the cell Peclet number exceeds the transition value. By the third assumption, the 
skew differencing formulae, Eqs. (A3. 4) through (A3. 6) are used but by the second 
assumption the flux contributions E3, N3 and B3 used in the calculation of the nine 
co-efficients are zero. By the first assumption, the flow direction switch at 
the node (I-1,J,K) is zero so that E2(I-1,J,K) is zero. As a result, the conditions 
in Eq. (A4.1) are satisfied and all nine coefficients are zero. Thus, the unknown 
values of the scalars at the unspecified opening are decoupled from the system of equa- 
tions . 

It is seen that the modifications to the flux contributions for the specified wall, 
axis of symmetry, and unspecified opening boundary condition are identical. 
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For a specified opening, the value of the scalar is known at the west face of the 
control volume for the node (l,J,K) since, as noted above, the physical boundary in 
this case is located at the west face of the scalar control volume. In general, the 
cell Peclet number is greater than the transition value whenever the boundary repre- 
sents a specified opening. Therefore, the flux contributions El, E2 and E3 at the 
node (l-l,J,K) are computed using the skew upwind differencing formulae, Eq. (A3. 4). 
However, since the boundary is coincident with the west face of the control volume for 
the node (l,J,K), there is no skewing of the flow at this face and, therefore, all E3 
at (I-1,J,K) are set to zero. 

The determination of the co-efficients for the node (l,J,K) in this case also re- 
quires the computation of north, south, back and front contributions. Since the speci- 
fied opening is coincident with the west face of the control volume for the node (I,J,K) 
and since this face is located midway between the nodes (I-1,J,K) and (l,J,K), it is 
necessary to double the value of the skewing interpolation factors used in calculating 
these flux contributions. Thus, for the specified opening: 

E3(I-1,J,K,L) = 0 L = 1,2, 3, 4 

(K 3 ) = min [ 1 , 2 (K 3 ) ] 

n n 

(K 3 ) = min [ 1, 2 (K 3 ) ] 

s s 

(K 3 )^ = min [ 1, 2 ( 1 ( 3 )^ ] 

(K 3 ) = min [ 1, 2 ( 1 ^ 3 )^ ] (A4.2) 

The analysis of the modifications to the flux contributions for a physical boun- 
dary at other faces of the control volume is similar and the results are summarized as 
follows . 
b) East Face 

(1) For a specified wall, axis of symmetry, or unspecified opening: 

E1(I,J,K) = 0 

E3(l , J,K,L) =0 L = 1,2, 3,4 

N3(I,J,K,L) = 0 

N3(I,J-1,K,L) = 0 
B3(I,J,K,L) = 0 
B3(I,J,K-1,L) = 0 
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(A4.3) 



(2) For a specified opening: 


E3(I,J,K,L) =0 L = 1,2, 3, 4 

(K 3 ) = min [ 1, 2 (K 3 ) n ] 

( K 3) s = m ^- n f 2(K3) g ] 

(K 3 ) fe = min [ 1, 2(K 3 ) b ] 

(K 3 ) = min [ 1, 2 ( 1 ( 3 ) £ 1 


(A4.4) 


c) South Face 

(1) For a specified wall, axis of symmetry, or unspecified opening: 


N2(I , J-l ,K) = 0 
N3(I,J-1,K,L) = 0 
E3(I,J,K,L) = 0 
E3(I-1,J,K,L) = 0 
B3(I,J,K,L) = 0 
B3( I , J ,K-1 ,L) = 0 

(2) For a specified opening: 


L = 1,2, 3, 4 


L = 1,3 


N3(I,J-1,K,L) = 0 L = 1,2, 3, 4 

(Ki) = min [ 1, 2(Ki) ] 

x e e 

(Ki ) = min [ 1, 2(Ki) ] 

w w 

(K]_) = min [ 1, 2(Kj_) b ] 

(Ki) f = min [ 1, 2(Ki) f ] 


(A4.5) 


(A4. 6) 


d) North Face 

(1) For a specified wall, axis of symmetry, or unspecified opening: 


N1(I,J,K) = 0 
N3(I,J,K,L) = 0 
E3(I,J,K,L) = 0 
E3(I-1,J,K,L) = 0 
B3(I,J,K,L) = 0 
B3(I,J,K-1,L) = 0 


L = 1,2, 3, 4 


L = 1,3 


(A4.7) 
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(2) For a specified opening: 


e) Front Face 
(1) For a 


(2) For a 


f) Back Face 
(1) For 


N3(I,J,K,L) =0 L = 1,2, 3, 4 

( Kl ) e = min [ 1, 2(Ki) e ] 

(K]_) = min [ 1, 2(Kx)^ ] 

(Ki) fe = min [ 1, 2(K 1 > b ] 

(K-^) f = min [ 1, 2(K 1 ) f ] 


(A4.8) 


specified wall, axis of symmetry, or unspecified opening: 


B2(l , J ,K-1) = ) 
B3(l,J,K-l,L) = 0. 
E3(I,J,K,L) = 0 
E3(I-1,J,K,L) = 0 
N3(I, J,K,L) = 0 
N3(I,J-1,K,L) = 0 


L = 1,2, 3, 4 

L = 1,2 

L = 1,3 


specified opening: 

B3(I,J,K-1,L) = 0 L = 1,2, 3, 4 

(Ko) = min [ 1, 2 (K 3 ) ] 

J e e 

(K 3 ) = min [ 1, 2 OC 3 ) ] 

J w w 

(Ki) = min [ 1, 2(Ki) ] 

n n 

(Kf) = min [ 1, 2(Ki) s ] 


(A4.9) 


(A4.10) 


specified wall, axis of symmetry, or unspecified opening: 


B1(I,J,K) = 0 
B3(I,J,K,L) = 0 
E3(I,J,K,L) - 0 
E3(I-1,J,K,L) = 0 
N3(l , J ,K,L) - 0 
N3(I,J-1,K,L) - 


L = 1 , 2 , 3 ,4 


L - 


1,2 
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0 


L - 1,3 


(A4.ll) 



(2) For a specified opening: 


B3(l , j ,K,L) = 0 
(K 3 ) = min [ 1, 

(Ko) = min [ 1 , 
J w 

(K-i ) = min [ 1 , 

A n 

(Ki)^ = min [ 1 , 


2(K 3 ) e ] 

2 ( K 3 ) u ] 

2<K 1>„ ] 
2( Kl ) s ] 


(A4.12) 


Ax_iaJ- JVel.o£i_ty 

The east and west faces of the axial velocity control volumes are co-incident with 
the vertical scalar gridlines. The axial storage locations of physical boundaries are 
the same as the axial locations of the axial velocities; the inverse is, of course, 
not true. The radial and lateral locations of the u velocity and scalar control 
volumes are identical. 
a) West Face 

Assume that a physical boundary is located to the west of the axial velocity control 
volume. If the physical boundary represents either a specified wall or an axis of sym- 
metry, then the axial velocity at this face is zero; if it represents a specified open- 
ing, then the axial velocity is a known, specified value. If the physical boundary 
represents an unspecified opening, then it is assumed that the axial velocity is known 
from the previous iteration. In every case, it is seen that the axial velocity is 
known at the west face of the axial velocity control volume. Therefore, the axial 
velocity at the west face of the control volume can be computed in the same manner as 
is used for any axial velocity interior node. 

The radial and lateral velocities for the west face of the axial velocity control 
volume are interpolated in the same manner as they are for the west face of a scalar 
control volume since the radial and lateral grid line locations for the scalars and 
axial velocity are identical. Thus, one can obtain the values of u, v and w and the 
flux contributions can be computed using Eqs. (A3.1) through (A3. 6 ) in the usual manner. 
Thus, it is seen, that the application of the boundary conditions requires no modifi- 
cation to the computation of the flux contributions. 
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b) East Face 

A similar analysis leads to the conclusion that no modifications are made to the 
computation of the flux contributions. 

c) South, North, Front and Back Faces 

Since the radial and lateral positions are identical for both the axial velocity 
and scalars, the modifications of the flux contributions for these faces of the axial 

velocity control volume are the same as those required for the corresponding face 
of the scalar control volume. 

Radijjl__V_eloci_t^ 

The radial velocity control volumes have radial and lateral positions identical to 
those for the scalar control volumes. 

a) South and North Faces 

No modifications are necessary. (See discussion for west and east faces of axial 
velocity control volumes.) 

b) West, East, Front and Back Faces 

Use modifications for corresponding face of the scalar control volume. 

Lateral Velocity 

The lateral velocity control volumes have axial and radial positions identical to 
those for the scalar control volumes. 

a) Front and Back Faces 

No modifications are necessary. (See discussion for west and east faces of 
axial velocity control volumes.) 

b) West, East, South and North Faces 

Use modifications for corresponding face of scalar control volumes. 

Special Treatment for Corner Regions 

At a corner, it is necessary to compute the flux contributions as a weighted 
average of the fluxes in the vicinity of both a wall and a fluid boundary (unspecified 
opening or specified opening). The weighting is equal to the fraction of the face of 
the control volume adjacent to the solid boundary. Since this procedure is not 
peculiar to the bounded Skew-upwind differencing procedure, but is part of the general 
procedures used for all finite-difference schemes in the TEACH computer program, it 
will not be described further. 
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A. 5 - The Bounding Scheme 


The calculation of the bounded skew-upwind differencing fluxes and, therefore, 
the determination of the coefficients to the finite-difference form of the equations 
of motion requires that the blending factor, y, be determined. The blending factor 
specifies the relative proportions of the flux computed using skew and upwind differen 
cing. For example, at the west face of the typical control volume, Eq. (A2.40) states 


Fw bsud YwFw SUD + ^ Yw ^ Fw UD 


(A5.1) 


The co-efficients including the local blending factors are listed in Table A-ll . The 
blending factors are limited to the range 0 < Y < 1. Examination of the terms in each 
of the co-efficients in Table A-ll lead to the following general conclusions: 

(1) For the principal co-efficients (those with a single subscript such as 
Ag, Ag, etc.), the term without a blending factor y is unconditionally 
non-negative and the five terms containing blending factors are uncon- 
ditionally non-positive. 

(2) For the edge co-efficients (those with two subscripts such as Ag E , A 
etc.), two terms are unconditionally non-negative and two terms are 
unconditionally non-positive. 

(3) The corner co-efficients (those with three subscripts such as.Aggg, 

Apsw> etc.) are unconditionally non-negative. 

Therefore, it is possible that some of the coefficients in Table A-ll will be negative 

It is desirable that all of the co-efficients of the finite-difference form of 
the equations of motion (Eq. (A1.7)) be non-negative for in this case the value of the 
dependent variable 4> at the node P is simply a weighted average of the values of 4 1 at 
the surrounding nodes exclusive of (the somewhat complicating) effects of local sources 
A bounding scheme is a procedure to limit the values of the co-efficients of the 
finite-difference equations in such a manner as to produce this physically realistic 
result. Its principal computational advantage is to exclude under- and overshoots of 
the solution during the iterative procedure; these oscillations can produce severe 
numerical instability. 
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The bounding procedure used herein is based upon the following sequence: 

(1) for each iteration, the solution for the distribution of the dependent variable 
<f> is obtained using the blending factors y determined during the previous iteration; 

the blending factors are- set to unity for the first iteration; 

(2) for each point P in the variable field, the maximum and minimum values for <t> 

are determined from the 26 neighbors: 

4> = max (<}>,) (A5 .2) 

max k 

(j: . = min (<{>,) (A5.3) 

mxn k 

where <j)^ is the value of <j) at any of the 26 neighboring nodes and k ^ P. 

It should be noted that the effects of sources are not included explicitly in calcu- 
lating <f> or <j> . since the explicit result is difficult to perform while avoiding 
max mxn 

double counting of the source effect. In any case, the results obtained will be no 
worse than to force the use of upwind differencing by making -y = 0. According to Pro- 
fessor A. D. Gosman, consultant to this project, the determination of 6 and 4> . 

max mxn 

is an area for which further research is needed. 

It should also be noted that one or more <j> may be excluded from the determination 

of <f> or (ji when the control volume under consideration is adjacent to certain types 
max mxn 

of physical boundaries. For example, if the x-axis is an axis of symmetry, then <j> 

u 

is excluded since it is identical to <£ and its use could yield an inappropriate range 

P 

of permissible values. 

(3) if d> . < d>_< d> , then the local value of the blending factor is unaltered 

mxn P max 

from its current value. However, if is outside of this range, then a new value of 
y must be determined. For this purpose each of the co-efficients in Table A-ll is 
written in the form: 

= A^ + yA^, k = N, S, NE, SE . . . (A5.4) 

where it has been assumed that the blending factors for each of the faces of the local 
control volume are equal. Then Eq. (A1.7) can be written as: 

^^[EA^ + y^A^] - Sp<j>p = E ^k^k + "Y ^ ^k^k + ^u (A5.5) 
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Now, if d>_ < 4> . , then from this equation, 
P mm 


Y > 


I Avd>i, — d> • Z A-/. + S + Sr,4' • 
KY k Y mxn k u P mm 


<t> . Z A k - I A-k^k 
mm K 


(A5.6) 


but if 4 > 

P 


max 


then 


Y 


Z A' 4, - 4 Z A,' + S + S 4 
k k max k u p max 


4 Z 

max 


" - Z A ?k 




(A5.7) 


Of course, in practice the value of Y is determined by use of the appropriate equality 
for Eq. (A5.6) or Eq. (A5.7). 

(4) In this fashion, the local value of -y(I,J,K) is determined. The co-efficients 
in Table A-ll are then recomputed using the blending factors and the new distribution 
of 4 is determined. For each face of the control volume at the point (I,J,K) there are 
two values of Y as determined for the two control volumes sharing the common face. 

Thus, a sufficient condition for assigning local values of the blending factors is 
given by: 


y = min [ y(I-1,J,K), y (I,J,K) ] 

Y e = min [ Y (I,J,K), y(l+l,3,K) ] 

Y = min [ y(l,J-l,K), y(I,J,K) ] 

s 

Y n = min [ Y (I,J,K) , y(l,J+l,K) ] 

Y f = min [ Y(I,J,K-1), y(I,J,K) ] 

y = min [ y(I,J,K), y(l,J,K+l) ] (A5.7) 

D 

This procedure for determining the blending factors effectively limits the value 
of the dependent variable range to values of its immediate neighbors. However, it is 
still possible that some of the main co-efficients will be negative. More restrictive 
schemes can be formulated that guarantee that all co-efficients of the finite-difference 
equations are non-negative but these procedures achieve this result by reducing the 
blending factor toward zero; the greater relative contribution of upwind differencing 
in this case produces larger amounts of numerical diffusion. Thus, the selection of 
a blending scheme represents a compromise between the undesirable effects of numerical 
instability and numerical diffusion. ^ 



The final blending factor field determines the ratio of the UDS solution 
present in the final solution. Since at present it is not possible to ensure 
that the blending factor field is unique for a given problem the solution 
generated by BSUDS2 is also not unique. Although it lies between UDS and SUDS. 
The reason for non-unique blending factor field is the following. 

The calculation of blending factors depends on the value of the coefficients. 
These coefficients are calculated on the basis of the flow field that existed 
in the previous iteration. Hence the initial blending factor field strongly 
depends on the first guess. If, due to a bad guess, for example, at some nodes 
in the first iterations some of the factors are set to zero, the following 
iteration will produce a solution for those nodes which will be in inbounds. 
Hence the blending factors will never be increased from zero subsequently, and 
the final solution will be nearer to UDS. On the contrary a good first guess 
will not cause these factors to be set to zero at the final solution would be 
nearer to SUDS. In short if a bad guess forces the introduction of more than 
the minimum required numerical diffusion in the solution, there is no 
mechanism for diffusion to be taken out in the final stages of the solution. 

In the present code this adverse effect is minimized by under-relaxing the 
blending factor. But this practice can some times cause a minor problem. In 
some cases it is possible to obtain convergence with the solution being 
unbounded at some nodes. This problem can be completely avoided by checking 
for the undershoots and overshoots after convergence has been obtained. If an 
unbounded solution is detected iterations can be continued until a bounded 
solution is obtained. The above approach to achieving a unique solution from 
BSUDS2 is only one of a number of options that can be used. Another approach 
would be to update the blending factor at every iteration and not under-relax 
them at all . 
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Nomenclature 


A 

a 

B1,B2,B3 

C 

D 

El ,E2 ,E3 
F 

1 
J 
K 

k l’ k 3 

L 

Nl ,N2 ,N3 

Pe 

r 

5 
s 
u 

V 

w 

X 

y 

2 
a 

r 

y 

Ax,Ay,Az 

6 

<f> 

P 
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Coefficient i n finite-difference equations (Eq. A1.7) 
Flow area 

Flux contributions, back face of control volume 
Convective flux 
Diffusive flux 

Flux contributions, east face of control volume 
Flux 

Index for nodes - axial direction 

Index for nodes - radial location 

Index for nodes - lateral direction 

Skewing interpolation factors 

Flux contribution index for B3, E3 or N3 

Flux contributions, north face of control volume 

Peclet number 

Radial coordinate 

Source 

Distance along streamline (e.g., see Fig. A-5) 

Axial velocity 

Radial (or vertical) velocity 

Lateral velocity 

Axial co-ordinate 

Vertical (or radial) co-ordinate 

Lateral co-ordinate 

Interpolation factor 

Exchange co-efficient 

Blending factor 

Internodal distances 

Azimuthal or co-ordinate 

Dependent variable 

Density 



U V w 

o , a ,0 


Subscripts 


B 

b 

BE,BN,BS,BW 

BNE,BNW,BSE,BSW 

BSUD 

CD 

E 

e 

F 

f 


FE,FN,FS,FW 

FNE,FNW,FSE,FSW 

k 


max 


min 

N 

n 

NE 

NW 

P 

S 

s 

SE 

sw 

SUD 

u 

UD 

W 

w 


Flow velocity direction switches 

Back node - see Fig. A-l 
Back face 

Back edge nodes - see Fig. A-l 
Back corner nodes - see Fig. A-l 
Bounded skew-upv?ind differencing 
Central differencing 
East node - see Fig. A-l 
East face 

Front node - see Fig. A-l 
Front face 

Front edge nodes - see Fig. A-l 

Front corner nodes - see Fig. A-l 

Dummy index 

Maximum value 

Minimum value 

North node - see Fig. A-l 

North face 

Northeast node - see Fig. A-l 
Northwest node - see Fig. A-l 

Node at center of control volume - see Fig. A-l 
South node - see Fig. A-l 
South face 

Southeast node - see Fig. A~1 
Southwest node - see Fig. A-l 
Skew upwind differencing 
Value not dependent on <f> 

P 

Upwind differencing 
West node - see Fig. A-l 
West face 
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7T 


f 

II 

Ilf 


Dummy node indicator 
Dependent variable 

Interpolated value 


Superscripts 

c ) u 
( ) v 
( ) w 
c )* 

( )' 

( ) m 

( ) , ( ) ,ef.c. 


Axial velocity 

Radial or vertical velocity 

Lateral velocity 

Transition 

For co-efficients, value not associated with y 
For coefficients, value associated with y 

Intermediate values of interpolated result for assumed velocity 
components 
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TABLE A-l 


Central Differencing Fluxes 


West 


w 


CD 


= [Pe w (l-a w ) + 1 ] <f> w + (a w Pe w "l) 4> p 


East 


J CD 


Dr 


P 0 (1-0 + 1 


+ ( a ePe e -D 4> E 


South 


S CD 


= [ p e s u-o + i 


+ ^sPeg" 1 ) <J> p 


North 


n CD 

D„ 


P P ^ d-otr.) + 1 !<}>„+ (a n Pe n ~D <i> N 


Front 


Pf 


CD 


Pe f (l-a f ) + 1 4 > f + (a f P ef -l) 4>j 


Back 


Pb 


CD 


* | P e b <1 -°b > * 1 1 + ( “b p e b -l) * B 
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TABLE A-2 


Upwind Difference Fluxes 


West 


^-p e L% * d-,“) ♦ 1 

D w e w l w W w J 


East 


e UD 


= P, 


u . 

0 „ ♦ » + <I>T 


South 


S UD _ [V v 


D c = Pe s +s + (1 -°» ) 


North 


n UD [ v , v, 

V ‘ Pe n °n *P * (1 -°„ ) *N 
N 


Front 


Ff 


UD 


= P e 


* F + ( 1 -f> * P 


Back 


' b UD f w . w, 1 

FT" Pe bl °b *P + (1 -°b > *b I 
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TABLE A-3 


Geometric Parameters Used in Interpolation Factors 







Nodes 

Used 

to 

Calculate : 


Face 

u 

V 

w 


Ax 



Ay 



Az 


East 

> 

0 

> 0 

> 0 

p 

to 

E 

P 

to 

S 

P 

to 

F 


> 

0 

1 0 

< 0 ' 

p 

to 

E 

P 

to 

S 

P 

to 

B 


> 

0 

< 0 

> 0 

p 

to 

E 

P 

to 

N 

P 

to 

F 


> 

0 

< 0 

< 0 

p 

to 

E 

P 

to 

N 

P 

to 

B 


< 

0 

> 0 

> 0 

p 

to 

E 

P 

to 

S 

P 

to 

F 


< 

0 

> 0 

< 0 

p 

to 

E 

P 

to 

S 

P 

to 

B 


< 

0 

< 0 

> 0 

p 

to 

E 

P 

to 

N 

P 

to 

F 


< 

0 

< 0 

< 0 

p 

to 

E 

P 

to 

N 

P 

to 

B 

North 

> 

0 

> 0 

> 0 

p 

to 

W 

P 

to 

N 

P 

to 

F 


> 

0 

> 0 

< 0 

p 

to 

W 

P 

to 

N 

P 

to 

B 


> 

0 

< 0 

I o 

p 

to 

W 

P 

to 

N 

P 

to 

F 


> 

0 

< 0 

< 0 

p 

to 

w 

P 

to 

N 

P 

to 

B 


< 

0 

> 0 

> 0 

p 

to 

E 

P 

to 

N 

P 

to 

F 


< 

0 

> 0 

< 0 

p 

to 

E 

P 

to 

N 

P 

to 

B 


< 

0 

< 0 

>_ 0 

p 

to 

E 

P 

to 

N 

P 

to 

F 


< 

0 

< 0 

> 0 

p 

to 

E 

P 

to 

N 

P 

to 

B 

Back 

2. 

0 

> 0 

> o 

p 

to 

W 

P 

to 

S 

P 

to 

B 


> 

0 

> 0 

< 0 

p 

to 

W 

P 

to 

S 

P 

to 

B 


> 

0 

< 0 

> o 

p 

to 

W 

P 

to 

N 

P 

to 

B 


> 

0 

< 0 

< 0 

p 

to 

W 

P 

to 

N 

P 

to 

B 


< 

0 

> 0 

> 0 

p 

to 

E 

P 

to 

S 

P 

to 

B 


< 

0 

> 0 

< 0 

p 

to 

E 

P 

to 

s 

P 

to 

B 


< 

0 

< 0 

> 0 

p 

to 

E 

P 

to 

N 

P 

to 

B 


< 

0 

< 0 

< 0 

p 

to 

E 

P 

to 

N 

P 

to 

B 
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TABLE A- 4 


INTERPOLATED VALUE FOR <J> - WEST FACE 

u w lO,V w > 0, all w w 

4> w m = ° w ^3 ^l^FSW + d~kl^ + (1 "°w ) ^3 ^l^BSW + ^ -k l^ 
+ ( 1 — k3 ) [k^^sw = (l - kl) 4>y] 


u w 1 °> v w < 0, all W w 

4 * T, , - k 3 [ k l‘ FKV , + * FM ] + ' l ~\? k 3 C k l*BNH * <1_k l ) ♦bu! 

* (1 - k 3> [XlV •' (l- k l> ♦„] 

U w < 0, v w > o, all W w 

Tm, = k 3 [k^ Fs + (1-ki) 4 F ] + U-o”) k 3 [k^ BS + (1-ki) ^ fi ] 

+ (l-k 3 ) [ki4> s + (1-ki) <j> p ] 

TJ W < 0, V w < 0, all W w 

O, = <\ 7 k 3 [ki<f>FN + (!-ki) 4> f 3 + (1-op k 3 [ki<t> BN + (l- k i) 4> B ] 

W W X W 

+ (l-k 3 ) [k 1 <|) N + (l-k^ <j)p] 


Then 


+ 

V 

++ 


V\ 

+- 

1 91 1 

- a 

<(> hi 

+ 

(l-o ) 

<J> m 

w ?n 

w 

w 


w 

W 

' ft 1 

V 

= a 

+ - 
i 

-&■ 

+ 

(l-a v ) 

in 

W 

w 

w 


w 

w 

w"' 

u 

= a 

C" 

+ 

(l-a u ) 

W 

*w ,n 



1 Ax Vw 
’ 2 Ay U w 

1 Ax W w 
’ 2 Az U w 
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TABLE A- 5 


INTERPOLATED VALUE FOR <j>-EAST FACE 


Ue - 0, Ve - 0, all We 


V’ = °e K 3 [ VfS + (1 - K 1 } V + (1 " a I ) K 3 [ VbS + V 


+ (1-K 3 ) [K^g + (1-K^ $ p ] 


Ue - 0, Ve - 0, all We 


V" = °e K 3 [K 1^FN + (1 - K 1 } V + (1 -< ) K 3 IVbN + (1 'V V 


+ (1-K 3 ) [K x ^ n + (l-K^tp] 


Ue < 0, Ve - 0, all We 


w 


V' ■ % K 3 [K 1*FSE + (1 - K 1 )4 FE ] + <K> K 3 'VbSE + <1_K 1 ) "bE 1 


+ a -V ‘Vse + (1 - k i h e) 


Ue_<_O x Ve_<_0_ L all We 
w 


■ "e K 3 [K 1*FNE + + “‘V K 3 [K 1 4 B K E + a 'V ^BE 1 


+ < 1 - K 3> [K lV + “‘W 


Then 


+ v ++ , v 4— 

d> tu = o <j> hi + (1—0 ) (j> hi 

e e e e e 

— v — t- v — 

<f> = 0 <t> II * + (1-0 ) cf) It' 

e e e e e 


u + ,, u N 

<j> ,,, = o <(> n • + (l-o ) 4> ti ' 
e e e e e 
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= min [1 


= min [1 


1 Ax 

2 Ay 

1 M 

2 Az 


Ve 

Ue 

We 

Ue 


|] 

II 
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TABLE A- 6 


INTERPOLATED VALUE FOR A-SOUTH FACE 


All Us . Vs — O, Ws - 0 


=JL 

C' = °S K 3 ‘VfSW + < 1 - R 1 >*SW 1 + (1 -°s ) K 3 1 VfSE + *SE ] 


+ (1-K 3 ) [K^pg + (1-K 1 )<t s ] 

Ail Us^_ V s _-_Oj_ Ws_<_0_ 

V- = °S K 3 l K !*BSW + (1 " K 1 H 'SW ] + (1 -°s> K 3 [ VbSE + *SE ] 

+ (1-K 3 ) [K 1 ^ bs + (1-K 1 )* S ] 

All Us_j_ V s _<_Oj_ W 3 _-_0_ 

V’ = °s K 3 [ VfW + (1_K 1 >+W 1 + (1_0 8 ) K 3 VfE + (1 - K 1 } *E ] 

+ (1-K 3 ) [K^ + (1-K 1 ) <^ p ] 

All P_ SjL Vs_<_0 JL Ws_<_0_ 

V- ■ % s ' K 1*BW + “'W + ( K> K 3 'VbE + (H 1> 4 E ! 

+ (1-K ) [K 1 * b + (1-KjHp] 

Then 

+ w ++ w H — 

' 0 s *,•" + (1 " 0 . )+ .- 
w — + w 

V" = °S V" + (1 "°s H s”'' 

V" = °S *8'" + (1 ’°s^s"' 

K ]_ - min l 1 , 2tt Ivt l ] 

K 3 = min [1> 2 ^ IvJ 11 



TABLE A- 7 


INTERPOLATED VALUE FOR (Jj-NORTH FACE 


All Un, Vn - 0, Wn ^ 0 


V" ■ x 3 'Vfv + + a -° U n> K 3 [K 1 »FE + (1 - K P *E ] 


+ <1-K 3 > [K 1 <t F + (1-Kj) ip] 


All Un. Vn - 0, Vn < 0 


V" ’ °n K 3 |K 1 *BW + V - “‘V K 3 [K 1 V + (1 'V +E 1 


+ <1-K 3 ) [Kj i B + (1-Kp) i p ] 

AU Un^ Vn_<_0^ Wn_-_0_ 

V" • °n K 3 lK l *FNW + (1 - K 1 ) W + (1 ”°n ) K 3 > K 1 *FNE + V 

+ (i-k 3 ) [k 3 t FN + <1-K ) ♦„] 


All Un, Vn < 0, Wn < 0 


V" ■ °l K 3 [K 1 <BNW + <1 - K 1 ) V’ + <1 -”n ) K 3 1K 1 W + V 


+ (1 - K 3> [K 1 + (1 - K 1 ) V 


Then 


+ 

w 

++ 

w 

<t- III 

= a 

<t> + 

(1-0 

n 

n 

n 

n 

_ 

w 

-+ 

, w 

1 ,ii 

= a 

$ , H + 

(l-o 

n 

n 

n 

n 


V 

+ 

, V 

,,i 

a 

+ 

-©■ 

(1-a 

n 

= n 

n 

n 

K = 

min 

[I, ^ 

1 Wn 

1 


1 ’ 2 Az 

' V n 

K = 

min 

[1. if 

,Un 

3 


2 Ax 

V n 


11 

|] 
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TABLE A-8 


INTERPOLATED VALUE FOR 4> -FRONT FACE 


lo_y_ all 

V” = °f <WfSW + (1 "V V + (1 - K 3 } [ VFS + (1 - K 1 } V } 

+ (l-op[K^ sw + (1-K^ ^ ] + (1-K 3 ) [K x 4> s + d-^) 4>pl } 
Uf_-_0^ V f _<0,_all_ w f 

(K 3 tK 1 + (1-Kj) ♦„] + Cl-V ! K ! *™ + <1 - K 1 > V 1 

- (l-0») (KjIKj ^ + (1-Ip ♦„) + d-K 3 ) [K x + (1-K a ) * p ]} 

Uf_<0»_Vi ~Pjl a 1 ! Wf 

V" * °£ lK 3 [K l "Vse + C1 - K l > V + <1 - K 3 ) 1K 1 *FS + (1 ‘V V’ 

+ (1-op (K 3 [K 1 * se + (1-Kj) * E 1 + (1 -KjMKj 0 S + (1-Kj) * p n 

U f _ < 0,_v f _ <0_^ all w f 

V. - »" « 3 ' K 1 * FN e + (1 - K 1 ) W + <1 - K 3 ) [K 1 *FK + <1 - K 1 > *F» 

+ ( 1 -op {K 3 [K 3 * ke + ( 1 -Kp * E 1 + < 1 -K 3 ) (K 3 + (1-Ip i p n 

Then 

+ v ++ -V. +- 

*f"' = °f *f'" + (1 • V *f ' " 

♦f.M = °f \ t u + (1 - 0 f } *f”' 

$f," = ♦j.n + (1 - °f) ^ f ... 

K, = min [1, 1/2 Azi V f M 

1 . >T7 U 

Ay Wf 

= min [1 9 1/2 AzjUf |j 
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TABLE A-9 


INTERPOLATED VALUE FOR 4-BACK FACE 


U b _ 2 _0^ »1£ “b_ 

V" - S 1K 3 i R l S« + (1 - K 1 > S 1 + IK T *S + <1 - K T ) S 11 

+ <K) (K 3 [K 4 bsu + d-K,) * BH ] + (1-K ) [K x 4 bs + (1-K 1 ) 4 b ]} 


S'" ■ S «3 !K 1 Sw + (1 - K 1> V + <1 - K 3 ) 1K 1 S + (1 - K 1> S !> 

+ W-S' {K 3 IK 1 *BNW + (1 -V Sw 1 + (1 - K 3> IK 1 Sn + (1 'V S !) 


V" ■ < (K 3 [K 1 +SE + (1 - K 1> S ! + <1 - K 3 ) [K 1 S + <1 - K 1 ) S 1} 

+ (1-0*) (K 3 [ Kl 4 ese + (1-Kj) * BE 1 + (1-K 3 ) [K x 4 bs + (1- Kl ) t. B )) 
D b ' °> v b " 0. *U w b 


S'" ‘ °b <K 3 lK l Se + (1 - K 1> S 1 + <1 - K 3 ) [K 1 S + <1 - K 1 ) S ]} 

+ ‘K* ,K 3 [k i She + a - K i> Se 1 + (1 ' K 3> [K i Sn + (1 'V S» 


+ v -H- ,, v s 4— 

S'" ■ S S'" + <1 ~S ) S'" 
S'" " °b V" + <1 '°b ) s’" 
S'" " °b S'" + <1 "% > S'" 


* 3 --” 
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TABLE A- 10 


SKEW UPWIND DIFFERENCING FLUXES 


West 


SUD * * * 

- = (Pe -Pep d, , + [Pe w (l-o^) + 1] % + (« w Pe w -1) 4 

W 


Pe w = 


: w 


A 

Pe = 


w u \ 
o - (1-a ) 

w w 


East 


Fe 


CTyn * . A * 

= (Pe -Pe ) 4 ,„ + [Pe (1-a ) + 1] 4 P + (« Pe “1) 4 r 

D eee e e r ee t 


Pe - -**■ , Pe = 

e C P e u / \ 

E o - (1-a ) 


South 


SUD * * 

- (Pe 8 -Pe s ) * s ,„ + IPe s U-o,) 


1] 4 C + (« Pe -1) 4 P 
S s s P 


D S * 1 

Pe - , Pe = 

s Cc s V . 

b o - (1-a ) 


North 


n 


SUD a a * 

— = (Pe -Pe ) 4 ,„ + [Pe (1-a ) + 1] 4 P + (a Pe -1) 4„ 
n. n n n’ n n P n n N 


N 


D N _ * 1 

Pe - ~r , Pe = 

n C-vr n v . _ * 

N a - (1-a ) 


n 


n 
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Frwt 


CTTT) & * * 

— — - (Pe f -Pe f ) + [Pe f (l-a f ) + 1] 4> F + (a f Pe f -1) tp p 

F 


D F _ * 1 

Pe r = 7~ •> Pe r ~ 

f Cp f w . 

a - (1-0 


Back 


qnn * •k * 

r 2 ■ v» + !P \ ( 1 -“b> + 11 4 p + ( “t v> 

B 


Pe = — , Pe 
b C B b 


% - (1 -“b> 
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TABLE A-ll 


COEFFICIENTS OF THE FINITE-DIFFERENCE EQUATION 

A (I,J,K) = B1(I,J,K) 

B 

- y ta U (1 -a W ) E3] (I,J,K,2) 

e e e 

- yAa-o") B3] (I, J,K,4) 

D b 

- Y [c v (l-a W ) N3] (I,J,K,3) 

n n n 

- Y t(l-a U )(l-o W ) E3] (1-1, J,K,2) 

wee 

- Y s [(l-a V )(l-a^) N3] (I-1,J-1,K,3) 

S n n 

A (I, J,K) = -y t (l-a U )(l-a W ) E3] (I,J,K,2) 

BE e e e 

- Y^ICl-o^d-o”) S3] (I,J,K,2) 

b b b 

- y [(l~a U ) a V (l-a W ) N3] (I,J,K,1) 

n n n n 

+ y [(l-o U ) (l-a V ) (1 -o W ) N3] (I,J-1,K,1) 
s n n n 

A (I, J,K) = -y [a U (l-o V ) (l-o W ) E3] (I,J,K,1) 

BN e e e e 

-Y b l(l-oJ)(l-oJ) B3] (I,J,K,3) 

-Y [(l-a V )(l-a W ) N3] (I,J,K,3) 
n n n 

+Y [(l-a U )(l-o V )(l-a W )E3] (I-1,J,K,1) 
w e e e 
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TABLE A-ll (Cont’d) 


COEFFICIENTS OF THE FINITE-DIFFERENCE EQUATION 

A(I,J,K) = - y t(l»o U ) (l-o V )(l-o W ) E3] (I, J,K,1) 

BNE e e e e 

- Y. [(1-a") (l-o^) (l-o”) B3] (I,J,K,1) 

b b b b 

-Y t(l"C U ) (l-o V )(l-o W ) N3] (I,J,K,1) 
n n n n 

A (I,J,K) = - y, [o. U (lV) (l-o") B3] (I, J ,K,1) 

BNw b b b b 

- y [o U (l-o V ) (l-a W ) N3] (I,J,K,1) 

n n n n 

+ y [a U (l-a V ) (l-o W ) E3] (I-1,J,K,1) 
w e e e 

A (I, J.K) = - Y [a U a V (l-o W ) E3](I,J,K,1) 

BS e e e e 

- B3] (I, J.K, 3) 

b b b 

+Y [(l-o U ) o V (l-o W ) E3] (1-1, J.K, 1) 
w e e e 

+ y [o V (l-a W ) N3] (I, J-1.K.3) 
s n n 

A BSE (I,J>K) = ~ Y e I(l-a“> a I ( 1 _a I ) E31 ^ I ‘ J » K ’ 1 > 

- yJ( 1"°“) a^(l-o^) B3] (I.J.K.l) 

b b b b 

+ y [(l-a U ) a V (l-a W ) N3] (I.J-l.K.l) 
s n n n 
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TABLE A- 11 (Cont'd) 


COEFFICIENTS OF THE FINITE-DIFFERENCE EQUATION 

a bsw ( i ' j ’ k) ' - V b t0 b °b <1 '°b ) B31 

+ y [o U 0 V (1-0 W ) E3] (1-1, J ,K,1) 
wee e 

+ y [o U a V (l-a W ) N3] (I,J-1,K,1) 
s n n n 

A ( I , J , K) = - y [a"(l-a”) B3] (I,J,K,2) 

Bw b b b 

- y to" a V (l-a W ) N3] (I,J,K,1) 

n n n n 

+ Y ta U (l-a W ) E3] (I-1,J,K,2) 

we e 

+ y [o U (l-a V ) (l-a W )N3] (I, J-1,K,1) 
s n n n 

A £ (I,J,K) = El (I , J ,K) 

- Y [(l-o U ) E3] (I, J,K,4) 

e e 

- Yja-O o” B3] (I , J ,K»2) 

b b b 

- y [ (l-a U ) a V N3] (I,J,K,2) 

n n n 

+ Y.[(l-o“) (l-a") B3] (I,J,K-1,2) 
r b b 

+ y [(l-o U ) dV) N3] (I, J-1,K,2) 
s n n 
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TABLE A-ll (Cont'd) 


COEFFICIENTS OF THE FINITE-DIFFERENCE EQUATION 

A (I, J,K) - B2(I,J,K-1) 

F 

u w 

- Y [o a E3] (I, J,K,2) 

e e e 

- y [a" a W N3] (I, J,K,3) 

n n n 

- y [(l-a U ) a W E3] (I-1,J,K,2) 

w e e 

+ Y f [\ B3] (I, J ,K-1,4) 

i b 

+ y U-a V ) o W N3] (I , J-1,K,3) 
s n n 

A ( I , J , K) - - Y t (l-a U ) a W E3] (I,J,K,2) 

FE e e e 

- y [ (1-° U ) a V a W N3] (I,J,K,1) 

n n n n 

+ yJ(1-0 B3] (I, J,K-1,2) 

r b d 

+ Y [ (l-a U ) (l-o V ) G W N3] (I,J-1,K,1) 
s n n n 

11 V u 

A (I, J,K) = - Y \o (l-o ) o E3] (I,J,K,1) 

FN e e e e 

- y [ d-a V ) a W N3] (I,J,K,3) 

n n n 

+ y [(l-o U ) (l-o V ) a W E3] (I-1,J,K,1) 
we e e 

+ B3 1 d,J,K-l,3) 

r d b 
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TABLE A- 11 (Cont'd) 


COEFFICIENTS OF THE FINITE-DIFFERENCE EQUATION 

A (I,J,K) = - Y t(l-o U ) (l-o V ) a W E3] (I,J,K,1) 

FNE e e e e 

- y [(l-a U ) (l-o W ) a W N3] (I,J,K.l) 

n n n n 

+ Y^t(l-o^) (l-o l) a* B3] (I j J ,K-1 , 1) 
f b b b 

A (I,K,K) = - y [o U (l-o V ) a W N3] (I,J,K,1) 

FNW n n n n 

+ y [a U (l-a V ) o W E3] (I-1,J,K,1) 
we e e 

+ Y r [°y (1-0^) B3] (I,J,K-1,1) 

r b b b 

A (I, J,K) = - y [a U a" o W E3] (I,J,K,1) 

FS e e e e 

+ y [(l-a U ) a V a W E3] (I,J,K-1,1) 
w e e e 

+ Y f o* B3] (I > J ,K-1,3) 

r b b 

+ Y [c V o W N3] (I , J-l ,K, 3) 

s n n 

A (I, J,K) = - Y [d-a U ) 0 V o W E3] (I,J,K,1) 

FSE e e e e 

- yJU-cO B3] (I,J,K— 1,1) 

f b b b 

- Y [ (l-o U ) a V a W N3] (I, J-1,K,1) 

s n n n 
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TABLE A-ll (Cont’d) 


COEFFICIENTS OF THE FINITE-DIFFERENCE EQUATION 
\sw (I ’ J ’ K) = V°e °l < E3] (I - 1 * J ’ K ’ 1) 

+ a l B3 ^ (i.J.K-i.D 

f D D D 

+ Y la U a V g W N3] (I,J-1,K,1] 

s n n n 

A (I, J,K) = - Y [o U o V o W N3] (I,J,K,1) 

FW n n n n 

+ y [o U o W E3] (1-1 , J,K, 2) 
wee 

+ Y f [o" or B3] (I,J,K-1,2) 

1 D D 

+ Y [o U (l-o V ) o W N3] (I,J-1,K,1) 
s n n n 

A (I,J,K) = N1 (I , J ,K) 

N 

- Y [o U (l-o V ) E3] (I , J,K, 3) 

e e e 

- Yv 1 ( 1-0 a Z B3] ( I » J » K » 3 > 

D D D 

- y [U-O N3] (I, J ,K,4) 

n d 

+ y [(1 -c U ) (l-o V ) E3] (I-l,J f K,3) 

wee 

+ y [(1-07) (1-0 B3] (I,J,K-1,3) 

I D D 
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TABLE A-ll (Cont'd) 


COEFFICIENTS OF THE FINITE-DIFFERENCE EQUATION 
A^CI.J.K) = - Y e [(l -a U e ) (lV) E3] (I,J,K,3) 

- Y k [(1-o^) (1-a^) a" B3] (I,J,K,1) 

b b b b 

- y [(l-a U ) (l-a V ) N3] (I,J,K,2) 

n n n 

+ Y r I(l-a“)(l-a7) d-a") B3] (I , J ,K-1 ,1) 
f b b b 

Vt ( i ’ j ’ k) = - V a b (1_a b ) a b B3] (I ’ J ’ Kjl) 

- y [a U (l-a V ) N3] (I,J,K,2) 

n n n 

- y l° U (l-o V ) E3] (1-1, J ,K,3) 

we e 

+ Y (1 -oT) a-o?) B3] (I, J ,K-1,1) 
f b b b 

A (I,J,K) « N2(I ,J-1,K) 

u 

- Y [o U o V E3] (I, J,K,3) 

e e e 

- a Z B3 ^ (I,J,K,3) 

b b b 

+ Y [(l~o U )o V E3] (1-1 s J,K, 3) 
w e e 

+ y.[al (lV) E3] (I,J,K-1,3) 
f b b 

+ Y to" N3] (I,J-1,K,4) 
s n 
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TABLE A- 11 (Cont’d) 


COEFFICIENTS OF THE FINITE-DIFFERENCE EQUATION 
A se (I,J,K) = - Y e [l-o") o V e E3] (I,J,K,3) 

- T b t (!"%)% °b B3] (I > J » K » 1) 

+ y [ (1— c!*) cl (l-o*) B3] (I,J,K-1,1> 
f b b d 

+ Y t(l-o U ) 0 V N3] (I» J-1»K,2) 
s n n 

A_.,(I,J,K) = - Y [o“ cZ o” B3] (I , J ,K , 1 ) 

SW b b b b 

+ y [° U o V E3] (1-1, J,K,3) 
wee 

+ Y Aol ol (l-o*) B3) (I » J ,K-1, 1) 

I b b b 

+ Y [a U o V N3] (I,J-1,K,2) 
s n n 

A ( I , J , K) = E2(I-1,J,K) 

W 

- yJ°Z a Z B3 3 (i.J.k.2) 

b b b 

- Y [o U a V N3] (I,J,K,2) 
n n n 

+ y [o U E3] (I-1,J,K,4) 
w e 

+ yAaZ (l-o*) B3] (I,J,K-1,2) 

r b b 

+ Y [a U (l-a V ) N3] (I,J-1,K,2) 
s n n 
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FIG. A-2 


GRID SYSTEM FOR SCALARS IN X-Y PLANE 
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FIG. A- 3 SCALAR CONTROL VOLUME IN X-Y PLANE 
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FIG. A-U STORAGE LOCATIONS FOR AXIAL (U) AND RADIAL (V) VELOCITIES 
RELATIVE TO SCALAR (0) GRID SYSTEM IN X-Y PLANE 



APPENDIX £& 




DESCRIPTION OF 3D-TEACH 


B.l Description of Aerothermal Model 

3D-TEACH is a computer code that can solve fully three-dimensional fluid 
dynamics problems. It can handle ax i symmetric, planar, or three-dimensional, 
elliptic, turbulent flows. It is one of a family of such codes with titles 
2D-TEACH, 2D(C)-TEACH, 3D-TEACH and, 2D-PREACH. 

The input to all of these codes is generalized such that different problems 
can be run without the need for Fortran programming between problems. Also, 
the physical models used can be turned on or off by input command. These 
collective features result in an extremely flexible system. 

The codes form a family in that: 

a) As computer codes they are written with the same format, menus and 
commands such that an operator trained to run 2D-TEACH can easily run 
3D-TEACH. 

b) All codes have an interactive nature using the IBM Conversational 
Monitor System (CMS), and use prompts and cautions to ensure smooth 
execution of a case. The operator has the choice of either CMS or 
Batch running. The recommended procedure is to set up a case and get 
it running on CMS, then switch it to Batch to complete execution. 

File modification and selection of running mode is identical for all 
codes. 

c) The same basic equations are solved and the same physical models are 
used, together with solution algorithms, in all codes such that 
regression is possible. This means that the same two-dimensional 
problem can be solved on 2D -TEACH and 3D-TEACH , and the same results 
will be achieved. The only difference between the 2D-TEACH and 
3D-TEACH is the additional dimension available in the 
three-dimensional code. Also, the post-processors available with 
3D-TEACH are necessarily more comprehensive than those with 2D-TEACH, 
in order to adequately display the results. 

The acronym TEACH (Teaching Elliptic Axisymmetric Characteristics 
Heuristically) represents a generic solution technique and these codes 
represent current production state-of-the-art calculations in terms of 
equations solved, physical models used, discretization of the equations, and 
solution algorithms. They are not perfect, but are a marked advance on 
one-dimensional flow calculations and global modeling that formed the previous 
capability. The structure of the codes has been made such that modular 
replacement can be carried out as better models and solution algorithms are 
developed, while the basic framework and operational features remain. 
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The 2D-PREACH code is similar in concept to the others, but uses different 
solution procedures. It is not considered further here. 

B.2 Outline of Calculation Procedure 

B.2.1 Equations to Be Solved 

The combustion process is an extremely complex turbulent flow. It is a 
somewhat daunting task to describe such a random flow of chemically active 
eddy structures in terms that can be currently solved and can provide useful 
answers for the designer of practical equipment. 

Currently, the most practical approach is to stay within the framework of 
continuum mechanics and to use a statistical description of the turbulence, 
coupled with the accepted Eulerian description provided by the Navier-Stokes 
equations of motion. Hence, an instantaneous quantity is described as the sum 
of a time-averaged value and a random, fluctuating value. 

When the statistical description of an instantaneous quantity is substituted 
into the Navier-Stokes equations (Reference B-l) and time averaged, the 
resulting equation set is known as the Reynolds equations (Reference B-2). 
These equations are similar to the Navier-Stokes equations except that 
time-averaged quantities are used, and for the appearance of time-averaged 
correlations of fluctuating quantities. 

Turbulent motions increase the apparent viscosity of a fluid by some orders of 
magnitude since there is a continuous transfer of energy from the mean flow 
into large eddies and thence, cascading down through progressively smaller 
eddies, to the molecular level where the energy is dissipated as heat. If 
laminar diffusion terms are therefore very much smaller than turbulent 
diffusion terms, then neglect of fluctuations in laminar viscosity is 
permissible. This results in simplification of the Reynolds equations. It is a 
frequently used practice (Reference B-3) to also neglect terms involving 
fluctuating density, although this implies that temperature differences in the 
flow are not large. This practice also results in simplification of the 
Reynolds equations. 

The simplified Reynolds equations are expressed in terms of time-mean 
quantities and also cross-correlations of fluctuating velocities such as 
uj'uj 1 . These terms are known as the Reynolds stresses, and result in a 
closure problem. Turbulence modeling provides the necessary descriptions of 
the Reynolds stresses in known or determinable quantities. When the flow 
consists of more than one chemical species, modeling is required also for the 
turbulent mass flux Uj'mi'. These terms arise from applying the 
statistical treatment of turbulence to an instantaneous species transport 
equation. The instantaneous energy equation is given the same treatment. 
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To model the turbulent mass fluxes it was assumed that, similar to molecular 
Schmidt and Prandtl numbers, there are turbulent Schmidt and Prandtl numbers 
that relate turbulent mass and heat diffusivities to momentum diffusivity. 
Closure to the Reynolds equations was provided by a particular turbulence 
model known as the two-equation or K- model. It relies on the eddy viscosity 
concept. Finally, the modeled equations were algebraically manipulated into a 
general form in cylindrical coordinates: 
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where: 


r eff ,0 


any of the independent variables 

an appropriate turbulent exchange coefficient, depending on 
what 0 represents 


= a so-called "source term" which lumps together all other 
terms in a given equation not included in the first four 
terms of Equation B.l . 


The equation given is for steady state flow. Reynolds averaging does result in 
the equations retaining time-dependent terms, but these have been dropped. 

This was done for two reasons: (1) compatibility with the present design 
system and (2) time averaging precludes dynamic behavior other than that 
induced deliberately through one of the independent variables. 

By way of example. Figure B-l gives the values of some of the items in 
Equation B-l. 
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Cg q is an additional constant* presently taken as being equal to C £l . 


Figure B-l Summary of Some of the Equations Solved in 3D-TEACH 
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B.2.2 Numerical Approach to Equation Solution 

The simultaneous set of main and auxiliary equations to be solved in a 
turbulent reacting flow with a liquid-fuel spray contains a significant number 
of individual equations, most of which are either ordinary or partial 
differential equations, and which are nonlinear. Numerical solution of these 
equations is necessary. Rearrangement into the general form represented by 
Equations B-l and B-2 enables one solution algorithm to be used for all 
equations. 

Conventional numerical methods available to solve equations of these types can 
be broadly divided into finite difference and finite element methods, although 
the dividing line is not distinct. Finite differences have a considerable 
background, and most solution approaches utilize this method. 

The finite difference analog of the differential equations is obtained by 
overlaying a computational mesh on the flow domain, and obtaining the basic 
finite difference form of the partial derivatives for every node of the mesh 
from a control volume approach (Reference B-4). The finite difference 
expressions, when substituted back into the differential equations, yield a 
set of linearized, algebraic equations for every node of the mesh. Thus, there 
are as many sets of equations as there are nodes in the calculation domain. 
These sets, along with the problem boundary conditions, can then be solved to 
give solutions for the entire flow field. 

Standard numerical techniques can be employed to solve the finite difference 
forms of the differential equations (Reference B-5). A steady-state implicit 
solution method is often used (Reference B-6); an initial guess is made of the 
field variables, and these guesses are iteratively updated until the solutions 
have converged. Convergence is deemed to have been obtained when the absolute 
sums of the residuals of each variable over the whole grid goes below a 
specified value. 

The relevant equations describing the flow motions, the physical models, and 
the solution techniques are assembled into the computer codes to carry out 
direct flow simulations on high-speed, large-core digital computers. Figure 
B-2 presents a flow diagram describing the calculation procedure, showing the 
assembly of the equations, the utilization of physical modeling, the computer 
solution, and the output for design use. The 3D-TEACH code conforms to this 
organization. 
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Figure B-2 Flow Diagram of Calculation Procedure 
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B.3 Solution Procedure 


With reference to Figure B-2, assembly of the equations governing the problem 
has been briefly described. The details of the computer solution of the 
resulting equation sets are to be described. 

Rearrangement of the equations into the general form represented by Equation 
B-l enables one solution algorithm to be used for all equations. The equation 
set is solved using a steady state, implicit, finite difference numerical 
procedure. Initial guesses are made of the field variables, and these guesses 
are iteratively updated until the solutions have converged. Convergence is 
deemed to have been obtained when the absolute sum of the residuals over the 
whole grid of each variable goes below a specified value. 

B.3.1 Discretization of the Equations 

The finite difference analog of the difference equations is obtained by 
overlaying a computational mesh on the flow domain to be calculated, and 
obtaining the basic finite difference form of the partial derivatives for 
every node of the mesh from a control volume approach, (Reference B-4). The 
finite difference expressions, when substituted back into the differential 
equations, yield a set of linearized, algebraic equations for every node of 
the mesh. To demonstrate, first in two-dimensions. Figure B-3 illustrates the 
mesh and the control volume established about a considered node, P. The 
control volume approach is based on the satisfaction of macroscopic physical 
laws such as conservation of mass, momentum and energy. The conservation 
property is essential when combustion is taking place. 
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Figure B-3 Control Volume for the Finite Difference Scheme 
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The code is written in both cylindrical and cartesian coordinates. The grid 
system consists of a set of coordinate lines intersecting in the x-r, x-©and 
r-Q planes for the cylindrical system. In the cartesian system the grid is 
formed by the intersection of x-y, x-z, and y-z plane lines. The intersections 
of these lines form the grid nodes at which all scalar properties are stored. 
Vector quantities are stored midway between nodes. Figure B-4 gives the finite 
difference grid control volumes for the scalar quantities and storage 
locations for the velocities in both coordinate systems. Note that compared to 
Figure B-3, there are two additional neighboring nodes, F and B, denoting 
Front and Back nodes in the z or© direction. The faces of the scalar control 
"volume are denoted by lower case letters. Figure B-5 shows typical scalar 
control volumes in perspective and gives the face areas and volume. Since the 
velocity components are located midway between the grid nodes, the control 
volumes for velocity components are formed by planes passing through the grid 
lines. Mote that since the control volumes for the velocity components are 
staggered (Figure B-6) the areas and volumes for these control volumes will be 
different from those of the scalar control volume. 

The finite difference form of the general partial differential equation is 
derived by supposing that each variable is enclosed in its own control volume, 
as illustrated in Figures B-3-6. The general 0 transport equation has a 
source term . This is expressed in linearized form and integrated over the 
control volume. The remainder of the transport equation is also integrated 
over the control volume, and added to the integrated source term. This yields, 

C £ <£ e ~ C W<*w + C H 0 n - c s 0 S + C B <b b - C F d> f » D E ( <* E - 0p) - D W ( 0p»4>w) 

+0 N ( <*fr V “ Qs ( <£ p - + °B ( “ *p> " °F ( - *F> * (Su + $p V 

B2. 


In the above equation the convection coefficients are defined as, 


C E * 5) e a e » c f s *)f a f 


and the diffusion coefficients are defined as: 


( LsHl±^) 

' A x / 



etc. 


a g ' a f etc. are the areas of the cell faces 
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Figure B-5 Perspective View of Scalar Control Volume, Giving Face Areas and 
Volume 
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Certain weighting factors are introduced into the variation of </> , the 
variable being calculated, and with the help of continuity. Equation B-2 can 
be manipulated and normalized to give the form, 

Aptf>p * A&s <t>H + A§ d>$ + Aj^W + Ap <*>p + Afi + S U 


where, B-3. 

Ap * Afj + As + A| + A^ ♦ Ap + Ag - Sp 

and Equation B-3 is the finite difference equation for 0 , and the main 
coefficients are defined as, 


An 

An 


Af| » 0 when P < »2 

a On - When -2 < P a 
2 e N 

88 Cn when P 

e N 


< 2 
> 2 


Similarly for A$. Ap etc., where the cell Peclet number is defined as, 

P * C N /0 M etc. 

®N 

There are several differencing schemes that can be used to evaluate the 
weighting factors. The values of the coefficients An, A$ etc. above were 
obtained using Spaldings Hybrid Differencing Scheme (HDS), of Reference B-7. 


The hybrid differencing scheme is unconditionally stable and the solution is 
bounded. It uses second order central differencing for convection and 
diffusion fluxes when the absolute value of cell Peclet number is less than or 
equal to two. When Peclet number is greater than two, first order upwind 
differencing is used for convection fluxes, and diffusion fluxes are neglected 
altogether. The switch of differencing is done both locally and directionally 
in the computational grid. Peclet number defines the relative importance of 
convective and diffusive transport. 

The finite difference Equation B-3, derived in the previous section could be 
used to obtain the velocity if the pressure field were known a priori. Since 
the pressure field is unknown, an iterative solution procedure, SIMPLE 
(Reference B-8) is used. SIMPLE is an acronym for J>emi Jjnplicit Method for 
Pressure Linked Equations. 

The essence of SIMPLE is that a pressure field is guessed, velocities are 
calculated from their finite difference equations, then the pressure and 
velocity fields are updated using a "pressure correction" equation which 
satisfies continuity. The procedure is repeated until the momentum and the 
continuity equations are adequately and simultaneously satisfied. The pressure 
correction equation can be derived from the continuity and momentum equations; 
the procedure is described below. 
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The finite difference form of the continuity equation can be written as: 


u) e a e - (p u) w a w + (P v) n a n - ( p v) s a s ♦ (pw) f a f - (pw) b a b = 0 ^ 
The momentum equations can be written as: 

A P Vp* - A n Vy* + A s v $ * + A e v e * + Ay Vy*- + Ap v F * + A B v b * + a$ (P s * - P p *) 
Ap Up* - Ay Uy* + A s u s * ♦ A e u e * + Ay Uy* + Ap up* + A b u b * + a w (P w * - P p *) 

A p w p* = a NWN* + A S W S* + a E w £ * + A w w w* + a FWF* + A B W B* + a b ( p F* “ p p*) 


In the above equations the pressure term has been separated from the source 
term and the (*) superscript denotes the values obtained from solving the 
momentum equations using the guessed pressure. An incorrect guess will give 
rise to a "mass source", Mp, in each cell because the continuity equation will 
not be satisfied. The mass source can be found by using Equation B-4. Hence 

Mp * ( p u*) e a € - (p u*) w a w + ( P v*) n a n - (p v*) s a s -(pw )paf + (pw ) b a b 


If the above equation is subtracted from Equation B-4., 


-Mp * ( p u' ) e a e - { P u‘ ) w a w + ( p v' ) n a n - ( P v* ) s a s -( pw ' )f a f + ( pw ‘ ) b3b 


where u e ! * (u - u*) e etc. 


B-5. 


The above velocity corrections can be calculated from the linearized momentum 
equations. 

Ap up* ■ a w (Py ! - Pp' ) 


B-6 
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Note that u p ' in the momentum equation control volume is u w ‘ for the 
continuity control volume and similarly for v p ' etc. On substitution of 
Equation B-6 in Equation B-5 and simplification, 

Ap Pp ' - A N P N ' + A S P S ' + A W P W ' + A E P£‘ + A f P f ’ + A g P g '+S u 

B-7 

where 


Aw ■ 

S u * - Mp 

Ap * A N + A$ + A w + Ag + Ap+Ae 

Equation B-7 is called the pressure correction equation which is solved to 
obtain corrected velocities and pressures. 


u e ■ u e * + Ug 1 
? e » P e * + P' 

etc. 

The difference equation for P' (Equation B-7) is in the same form as the 
difference equations for 0 (Equation B-3) an d ence a single solution 
algorithm can be used to solve all difference equations embodied in the 
numerical method. 

Since the SIMPLE procedure computes the variable fields successively it is 
highly flexible with respect to the methods of solution which it will admit 
for the difference equations. At present the following line by line iteration 
method is employed. This method is also known as Alternating Direction 
implicit Method (Ref. B-9). The ADI methods were Initially formulated for 
unsteady equations; their adaptation to steady state equations is sometimes 
also known as Alternating Direction Herative Methods. 

The finite difference Equation (B-3) to be solved is 


Ap *p » A N * N + A$ 0$ + A W + A E +A F 0 F + A B 0 B +S U 
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where 0 stands for u, v, p, K, e , and H successively. This equation can be 
recast in the following form 


Ap 0p * Af^ 0 N + *5 0 S + e ' 
or 

Ap 0j - An <f>j+l + As 0j-l + C j* 


B-8. 


To solve the equations for points on each line (e.g., N-S line) values on 
neighboring lines are assumed to be temporarily known. The equation for each 
point on the N-S line then reduces to one where only three values ( 0p, 

0 N> 0 S 1n Equation B-8) are unknown. An equation of this type can then be 
solved by the Tri-Diagonal Matrix Algorithm (TDMA), which is explained below. 

Equation B-8 can be rearranged for the jth point as 


0j - Bj 0j + l + Cj 0j-l + °j 

where 

Bj - A N /Ap , Cj * A$/Ap 

Dj * (Ay « w + A£ 0£ + Ap0p + ^ S u>/Ap 


The points on the computation grid range from 1 to Nj in the N-S direction 
with points 1 and Nj on the boundaries. Since the boundary values <b \ and 
0Nj are known, equations for 02 to 0 NJ-l are solved. The set of 
equations then becomes: 


*2 * B 2 03 + C 2 0 1 + d 2 
0 3 * B 3 0 y + C 3 0 2 + o 3 


0 Nj-l 


B Nj-l 



+ C 


Nj - 1 0 Nj -2 * 


°Nj-l 


B -9 
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Now, since 01 is known, 02 can be eliminated from Equation B-9 and so on, 
yielding a general recurrence relation 


8 Aj *J+1 + D j‘ 

B-10 


To get the relation for Aj and Dj' Equation B-10 is written as 
♦ j -1 “ Aj-l * j ♦ D j -1 

Now putting in the value of 0 j 


a N 


\A P - A S . Aj_i 



+ 


Cj l * A $ Oj.i 
. A p - A S 


B— 1 1 


Comparing Equation B-10 and B-ll yields coefficients for the recurrence 
formula 


where A j * (V < A P ' A S A j-1 > ) j 


B-12 


D j ' (f A S D j-1 + c j ) / {Ap - A S Aj _i ) ) . 

B- 13 


Using Equations B-12 and B-13, 0 j can be calculated from Equation B-10. 
Having solved for 0 j on one N-S 0 j> s on the next N-S line are solved 
and so on until the entire solution domain is swept. The same treatment is 
then applied in the W-E direction and finally, in the F-B direction. It is 
usually necessary to sweep between 1 and 3 times per iteration for optimum 
solution time. 

The coefficient matrix formed by the finite difference equation of 0 should 
satisfy the stability condition. 


Now 


Ap i* I A„ I 
A P * ^ A n - SP 
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So if 


S 0 < Q 


the above criteria is satisfied. In the solution procedure care is taken so 
that Sp is always less than or equal to zero. 

In the process of the computations, convergence is assessed at the end of each 
iteration on the basis of the "Residual Source" criterion. The residual source 
R 0 is defined as 

*♦ * *p A „ % ' S 

It is required that: 

2 I I < « R$ R e f 

for each finite difference equation. 

R 0 Ref is the fixed flux of the relevant extensive property fed into the 
domain of calculation, and e is of the order of 10 - 3 . 

When it is the equations for mass fraction of species that are being solved an 
additional convergence criterion requires that the sum of the mass fractions 
at each node is ^ (1 +e). 

When the flow is of variable density it is initially required that the change 
in density in one iteration at every node must also be less than £ 

or 



Since the enthalpy values in the calculation domain do not conform with the 
species mass fractions during the first few iterations, temperature and 
density are not updated for the first 10 - 25 iterations. If the density 
gradients are steep, density is updated every second or third iteration after 
the first update. 

A typical convergence plot is shown in Figure B-7. 
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Figure B-7 Typical Convergence Plot 


Since the finite difference equations are nonlinear in nature the convergence 
is facilitated and sometimes divergence is avoid e by under-relaxing the value 
of being calculated as: 


R , + { 1-F) <Pp 01 d 

<t> p y ' 

B-14 

where F is an under-relaxation factor which is less than one. 

The way in which the above relation is introduced into the numerical procedure 
is as follows: 

*p * Ap/F 

Sj ■ ♦ o - n Ag pOid 

B— 15 
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It can easily be shown that the effect of introducing the above modifications 
is to under-relax 0 p according to Equation B-14. From Equation B-3 we have 


0 

P 



* $u 


Ap 


B-16 


putting in the under- relaxation factors, 


4> 


R 

P 


£ A n *n + S U R 

A* 

P 


B- 17 


putting the value of A p R and S U R from Equations B-14 and B-15 in 
B-16, 


I A n ♦ „ * S u ♦ (l-F) A* ^011 
K . n P P 

*P ' A* 

P 

R . (n A " *" * M F * H-F) Old 
*p \ — ^ 

, R » <p Mew f + (i-F) 0 01 d 
P P P 


B-18 


B-19 


B-19 


It should also be noted that the effect of under-relaxation is to make the 
coefficient matrix more diagonally dominant. 
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The various steps in the numerical procedure can now be summarized as follows: 

1. Guess fields for all variables. 

2. Assemble coefficients of momentum equations and solve for U* and V* 
using prevailing pressures. 

3. Solve the pressure correction equation and update velocities and 
pressures. 

4. Solve equations for other variables. 

5. Update fluid properties such as viscosity and density. 

6. Test for convergence. If not attained use prevailing fields as new 
guesses and repeat from step 2 until convergence is attained. 

In general, it is necessary to specify 0 or its gradient at the boundaries 
of the calculation domain. There are six types of boundaries that can be 
assigned to a fluid block: 

1. Plane/axis of symmetry 

2. Unspecified wall 

3. Specified wall 

4. Unspecified opening 

5. Specified opening 

6. Periodic 

A specified boundary is one for which all boundary values such as velocities, 
temperatures, etc. are given. An unspecified boundary is one for which the 
boundary values are calculated by the code. If the flow structure is 
repeating, periodic can be used to compute only one segment. 

On a plane/axis of symmetry the gradient of all 0's except v-velocity is put 
to zero; v-velocity itself is set to zero. It should be noted that the 
symmetry condition should be applied to all block faces which constitute the 
plane/axis of symmetry. Most walls are specified, with all velocities set to 
zero (no slip condition). A moving wall is modeled with the no slip condition 
by specifying nonzero velocities in the plane of the wall. A porous wall is 
modeled by specifying nonzero velocities normal to the plane of the wall. At 
the outflow, it is required that there be no negative axial velocity 
components. At high Reynolds numbers this requirement makes specification of 
boundary values of all 0 except u redundant. The axial velocity is specified 
thus, 


u n l, j * u ni-l , j + U I NC 


B-21 


where ni is the outflow boundary and Utmq is calculated such that the total 
mass outflow is equal to mass inflow. Alternatively, if an exit velocity 
profile is known, it can be specified. 
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The calculation mesh is constrained by the coordinate system, which presently 
has been selected as orthogonal. Therefore, curvilinear geometries have to be 
represented in the form of discrete steps or "staircases." The specified 
blockage boundary condition permits this representation and allows inflow and 
outflow through elements of these staircases. In addition, the condition 
allows solid bodies to be placed inside the flowfield and to contain mass 
sources or sinks within them. This capability is written in generalized form 
and confers considerable geometric flexibility on the code without the need 
for interproblem reprogramming. 

Adjacent to solid boundaries, the local Reynolds number of the flow based on 
local velocity and distance from the wall becomes very small and the 
two-equation turbulence model, which was developed for high Reynolds numbers, 
becomes inadequate. Although a version of the two equation model that can 
handle both high and very low local Reynolds numbers exists (Reference B-10), 
its application requires a large number of grid nodes (more than 30) in the 
wall layer. This is due to the steep gradients of properties in the wall 
region (Reference B-ll). 

To avoid these difficulties, it was argued that the flowfield in the 
calculation domain is not influenced to first order by the details of the flow 
at the walls (Reference B-12). Consequently, as a matter of computational 
efficiency and economy, the high Reynolds number version of the turbulence 
model was retained and a Couette-flow analysis was used to give an equilibrium 
boundary layer on all solid surfaces bounding the calculation domain. The 
resulting wall functions are used to link the walls to the near-wall nodes of 
the finite difference grid. 
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